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In this paper we theoretically discuss how quantum simulators based on trapped cold bosons 
in optical lattices can explore the grand-canonical phase diagram of homogeneous lattice boson 
models, via control of the trapping potential independently of all other experimental parameters 
(trap squeezing). Based on quantum Monte Carlo, we establish the general scaling relation linking 
the global chemical potential to the Hamiltonian parameters for the Bose-Hubbard model in a 
parabolic trap, describing cold bosons in optical lattices; we find that this scaling relation is well 
captured by a modified Thomas-Fermi scaling behavior - corrected for quantum fluctuations - in the 
case of high enough density and/or weak enough interactions, and by a mean- field Gutzwiller Ansatz 
over a much larger parameter range. The above scaling relation allows to control experimentally the 
chemical potential, independently of all other Hamiltonian parameters, via trap squeezing; given 
that the global chemical potential coincides with the local chemical potential in the trap center, 
measurements of the central density as a function of the chemical potential gives access to the 
information on the bulk compressibility of the Bose-Hubbard model. Supplemented with time-of- 
flight measurements of the coherence properties, the measurement of compressibility enables one to 
discern among the various possible phases realized by bosons in an optical lattice with or without 
external (periodic or random) potentials - e.g. superfiuid, Mott insulator, band insulator, and 
Bose glass. We theoretically demonstrate the trap-squeezing investigation of the above phases in 
the case of bosons in a one-dimensional optical lattice, and in a one-dimensional incommensurate 
superlattice. 
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I. INTRODUCTION 

The impressive recent advances in the engineering of 
interacting Hamiltonians for cold trapped atoms suggest 
the possibility of experimentally determining the equi- 
librium (and out-of-equilibrium) behavior of fundamen- 
tal theoretical quantum many-body models. In partic- 
ular experiments on cold atoms in optical lattices have 
demonstrated the ability of implementing fundamental 
lattice models (of Hubbard type) with full tunability of 
all Hamiltonian parameters, as well as of the lattice ge- 
ometry and dimensionality [U |5j . Remarkably, the un- 
derstanding of a large class of lattice many-body models 
(such as interacting fermions, frustrated quantum mag- 
nets, etc.) keeps defying standard theoretical and compu- 
tational approaches. Hence, the perspective of realizing 
analog quantum simulators [3] based on cold atoms, liter- 
ally implementing the physics of the challenging models 
in question, represents a most promising and innovative 
route towards their understanding. 

Nonetheless, several obstacles still separate the cur- 
rent experiments from giving original answers to long- 
standing questions in quantum many-body physics on a 
lattice. One problematic feature of cold-atom systems in 
a lattice is thermometry [511151 1161 .'' ] : given that trapped 
cold atoms are not coupled to a thermal reservoir, their 
entropy can be controlled but not their temperature, and 
the knowledge of the temperature as a function of en- 
tropy is a problem which requires the prior knowledge 
of the equation of state of the many-body system under 
investigation [TH [T71 [TH] ■ A second puzzling feature of 



current cold-atom experiments - which will be the main 
concern of this paper - is the intrinsic inhomogeneity 
of these systems, imposed by the presence of an over- 
all parabolic trapping potential. While the size of the 
lattices realized in experiments can easily beat the capa- 
bility of current classical simulation methods (at least in 
dimensions d > 1), the parabolic potential significantly 
limits the size of the lattice over which a uniform phase 
of the Hamiltonian model is realized. Indeed, a parabolic 
trapping potential V{i) = Vt{ri — rg)^ (where is the 
position of the i-th lattice site and tq is the center of the 
trap) imposes a site dependent chemical potential and 
site-dependent energy jumps AE{i) = V{i + 1) — V{i) = 
Vt{ri+i — ri){ri+i + ri — 2ro) between nearest-neighboring 
sites. If the local many-body phase at a given filling 
n, realized around i-th site at the local chemical poten- 
tial — V{i), is not protected by a particle/hole gap 
Ag ^ AE{i), its existence is necessarily restricted to 
a narrow neighborhood of site i. This aspect introduces 
significant limitations in the capability of realizing phases 
which are particularly sensitive to the filling, and which 
display a small or no particle/hole gap over the ground 
state (see e.g. the supersolid phase of strongly correlated 
bosons [SHZ] or the Bose-glass [5] for bosons in a random 
potential) . 

In principle one could consider the coexistence of sev- 
eral different phases in the trap as an advantage, given 
that one single experiment samples a significant portion 
of the phase diagram of the Hamiltonian implemented 
in the system via the chemical potential modulation im- 
posed by the trap. [121 120] Unfortunately this form of 
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"parallelism" is generally hard to enjoy, as most mea- 
surement protocols for trapped cold atoms employed to 
date give access to global observables, which collect a 
signal from all the different regions of the trap. This 
fact becomes particularly inconvenient when considering 
the measurement of the excitation spectrum, as sampled 
e.g. via lattice-modulation and two-photon Bragg spec- 
troscopy [21 24 . In fact it is quite difficult to associate 
the different contributions to the excitation spectrum 
with different spatial regions of the trap, a fact which 
makes it very hard to extract precise information on the 
structure of the excitations associated with a particular 
phase of interest. In particular, the presence of a trap 
imposes the existence of a halo of dilute particles at the 
cloud boundary, which in principle can always host low- 
energy excitations. Consequently, the fundamental ques- 
tion of the presence or absence of an energy gap over the 
ground state of a particular phase, realized locally in the 
trap, becomes a formidable task for current cold-atom 
simulators. 

It is important to mention that single-site address- 
ability in optical lattices is becoming reality in recent 
experiments on quantum gas microscopy p5H2^ . This 
achievement definitely allows one to enjoy the above men- 
tioned parallelism in the chemical potential exhibited by 
trapped experiments |19j . but only at the level of local 
properties, defined on spatial regions over which the vari- 
ation of the local chemical potential can be considered 
weak; and only as long as the so-called local-density ap- 
proximation can be fully trusted. Yet non-local proper- 
ties of homogeneous phases, including correlation func- 
tions and collective excitations, cannot be extracted from 
an inhomogeneous sample, except for those associated 
with the region of maximal homogeneity, namely the trap 
center. 

In this paper we propose an experimental protocol 
which aims at circumventing most of the difficulties listed 
above, while taking advantage of the trap as a low-energy 
probe for the properties of the system (trap-squeezing 
spectroscopy [30]). The fundamental idea relies on the 
fact that trap effects are minimal in the center, both 
V{i) and AE{i) vanish, so that a local, homogeneous 
phase can be established over a significant portion of the 
system, roughly of the order of ^ 10 lattice sites in each 
spatial direction. This sizable portion of the lattice can 
be now regarded as the system of interest, while the rest 
of the lattice as the environment, acting in particular as a 
particle reservoir (see Fig. [l] for a sketch) . In light of the 
above discussion, we can conclude that the trap center 
realizes the "textbook" quantum simulation of the par- 
ticular Hamiltonian implemented in the system in the 
grand canonical ensemble (namely at a nearly uniform 
local chemical potential). What remains to be shown 
is how to control the crucial parameter of the quantum 
simulation at the trap center, namely its local chemical 
potential, and how to retrieve selective information on 
the local collective phase realized there. 

We here specialize our discussion to the case of bosons 
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FIG. 1: An atomic cloud in a smoothly varying trapping po- 
tential - such as a parabolic trap - can be virtually divided 
into two regions: a central C region over which the density 
is essentially uniform (within a tolerance e), and whose lo- 
cal chemical potential I-l{N, Vt) essentially corresponds to the 
global chemical potential of the atomic cloud; and a periphery 
region C, in which the density varies rapidly in space, due to 
the trapping potential. The nearly homogeneous central re- 
gion can be regarded as the "system" in a grand-canonical 
ensemble, while the periphery can be regarded as a particle 
"reservoir". Increasing the trapping potential (as indicated 
by the black arrows) increases the chemical potential of the 
reservoir C, and, at equilibrium, that of the system C. Hence 
the chemical potential of the system can be controlled contin- 
uously by the trapping potential Vt. 



in optical lattices, governed by the Bose- Hubbard Hamil- 
tonian in an external potential, and leaving the case of 
fermions to future work. In particular we quantitatively 
discuss how, tuning the trap frequency independently 
of all others experimental parameters, one gains direct 
access to the control of the chemical potential in the 
trap center. Remarkably, numerical simulations on the 
Bose-Hubbard model on the hypercubic lattice show that 
there exists a simple relation between the strength of the 
trapping potential and the chemical potential: this rela- 
tion follows approximately the predictions which can be 
obtained both from the atomic limit in a lattice, and 
from the Thomas-Fermi theory for weakly interacting 
gases in continuum space, with deviations due to quan- 
tum and lattice corrections. Remarkably, such deviations 
are quantitatively captured at the mean-field level. This 
means that, unlike the case of the temperature \Tl\ 118) . 
the accurate knowledge of the chemical potential of a 
strongly correlated bosonic system in an optical lattice 
might not require in general an extensive ab-initio calcu- 
lation. This makes the trapping potential a fundamental 
experimental knob for the quantum simulation, whose ef- 
fect on the system's parameters can be readily assessed. 

Nonetheless, the tuning of such a knob has to be done 
very carefully, and in general it cannot be done after 
loading the atoms in the optical lattices. In fact, increas- 
ing the strength of the trapping potential at equilibrium 
leads to the transfer of particles from the wings to the 
center of the trap, but the tunneling amplitude associated 
with such a transfer can be extremely small in strongly in- 
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teracting systems, requiring then exceedingly slow ramps 
of the trapping potential to enforce adiabaticity in the 
protocol. A simple way to circumvent this problem is 
proposed, based on the loading of atoms in the trap in 
the weakly interacting regime (namely without optical 
lattice) followed by the ramp of the optical lattice. 

Once full control on the chemical potential at the trap 
center has been achieved, selective information can be 
retrieved on this region of space by microscopy of the 
atomic cloud. We here propose a measurement scheme of 
the average central density based on two tightly focused 
crossed beams resonant with different optical transitions 
of the atoms; as already mentioned, other quantum gas 
microscopy schemes (with resolution as high as one lat- 
tice spacing) have been proposed or even become exper- 
imentally available in the recent past. The knowledge of 
the central density as a function of the chemical potential 
provides then the fundamental information on the bulk 
compressibility of the Hamiltonian model implemented 
in the system, and hence on the particle/hole gap over 
the ground state. We present the application of trap- 
squeezing spectroscopy to the measurement of the bulk 
phase diagram of the Bose-Hubbard model in d = 1: in 
particular we discuss the case without any applied exter- 
nal potential, and the case in which an incommensurate 
potential with two wavelength components is applied to 
the system, giving rise to an extended Bose-glass phase 
whose compressible nature is perfectly captured via trap- 
squeezing spectroscopy. 

The structure of the paper is then as follows. Section 
[m introduces the model investigated, the local-density 
approximation and the fundamental role played by the 
average central density in the trap; Section [HI| discusses 
the relationship between the chemical potential in the 
trap center and the trap strength for the Bose-Hubbard 
model in dimension d = I, 2, and 3, with or without 
an external potential; Section |IV| discusses the issue of 
equilibrium preparation of the system at a given trap 
strength, based on a two-step protocol for the ramp of the 
optical potentials; Section |V] discusses a proposal for the 
selective measurement of the average density in the trap 
center; applications to the Bose-Hubbard model without 
and with external superlattice potentials are presented in 



Section VI and |VII| respectively; Section VHI compares 
the central compressibility with the global compressibil- 
ity of the atomic cloud, recently measured in experiments 
[ST] : and finally Section IX is devoted to conclusions. 



II. LOCAL DENSITY APPROXIMATION, 
AVERAGE CENTRAL DENSITY AND CENTRAL 
CHEMICAL POTENTIAL 

In this section we focus our attention to the general 
case of the Bose-Hubbard model on the d-dimensional 
hypercubic lattice in an external potential, composed of 
a rapidly varying part, given by a superlattice potential 
created by an additional standing wave applied to the 



system [351 133] : a-nd a slowly varying parabolic part, 
coming from the overall gaussian profile of the lasers ap- 
plied to the system. 



niJ, U, V2, Vt) = HoiJ, U, V2) + 14 ^(r 



(1) 



noiJ,U,V2)^ - h.c.) +^^n,(n,-l 



where 



gi{{ai}, {4'i}) = ^ cos^(27rai ii + (pi) 
1=1 



(2) 



(3) 



is a one-color superlattice potential. Here ii {I — 1, d) 
are the coordinates of the i-th site and (ij) are the pairs of 
nearest neighbors on the d-dimensional hypercubic lattice 
of size L'^. 

Experiments on cold atoms in optical lattices are typ- 
ically performed with a fixed number of particles N. 
Nonetheless, given that the particle number is a good 
quantum number of the Hamiltonian, at T = the sys- 
tem will be in a definite A'^-sector even in the grand- 
canonical ensemble. This allows us to regard the canon- 
ical system with N particles as equivalent to a system in 
the grand-canonical ensemble with Hamiltonian 



— H 



t^ 



(4) 



where /i — ^{Vt, N, J,U,V2) is the chemical potential 
which establishes N particles in the ground state of the 
Hamiltonian 'H{J,U,V2,Vt). Moreover we can rewrite 
Eq. Q as 



"H^ = Ho -^n{i)r 



(5) 



where we have introduced the site-dependent chemical 
potential fi{i) = fi — Vt{ri — ro)^. If /i(i) is a slow- varying 
function in space around the reference site i* over the 
typical length scale given by the correlation length of the 
Hamiltonian H^^i*) — Ho — Ai(**) X^i "■jj '^^^ adopt 
the local-density approximation (LDA): this amounts to 
considering that the trapped system behaves around the 
site i* in the same way as its bulk conterpart at the homo- 
geneous chemical potential fi = ^i{i*). In particular, in 
absence of a superlattice potential, the local density (n^) 
for i Ki i* will be extremely close to the homogeneous den- 
sity in the ground state oiH^[i-'y In presence of a super- 
lattice potential, the local density (n^) will be very close 
to the density of the bulk system at a point experienc- 
ing the same superlattice potential; moreover the average 
density over a period of the superlattice (or quasi-period 
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for incommensurate superlattices) centered around the 
site i* will approximate very well the average density in 
the ground state of the bulk system. All these expec- 
tations for the behavior of the density are fully verified 
by numerically exact calculations on the Bose-Hubbard 
model with or without a superlattice [34l - [36j . The LDA 
typically breaks down when the local chemical potential 
fi{i) approaches a critical value sitting at the boundary 
between two phases in the bulk system, so that the cor- 
relation length of the bulk system diverges. 

The validity of the local-density approximation for 
parabolic traps implies that the trapped system faithfully 
probes the density of the bulk system at many different 
values of the chemical potential away from critical points. 
As mentioned in the introduction, this form of "parallel" 
sampling of the bulk phase diagram is only valid for lo- 
cal properties, and it cannot be exploited experimentally 
unless single-site addressability is achieved. On the con- 
trary, imaging of the atomic cloud over a length scale R 
of 5 — lO/Ltm (corresponding to ^ 10 — 20 lattice sites of 
an optical lattice with A ~ 800 nm) can be achieved more 
conventionally via large-aperture optics |37) (see Section 
|v]for a more detailed discussion). This means that the 
information on the local density can be retrieved if (rii) 
does not change appreciably over the lengthscale R. For 
a weak enough trapping potential, such a condition can 
be easily met at the trap center, where the variation in 
the local chemical potential is the slowest. We hence 
introduce the average central density 

nc =■ ^ Yj^n.) (6) 

where |C| is the size of the C region built around the 
trap center zq- The C region is then defined as verifying 
a condition of quasi-homogeneity 

k^czi!^ < e (7) 

with e ^ 1. Another important observation singles out 
the trap center as the most interesting region of the trap. 
Indeed, even when it is possible to measure the density 
at the single-site level, one can exclusively achieve the 
knowledge of a local observable Ai which, via the LDA, 
can be associated with that {A) of the bulk system as a 
function of the chemical potential, A{p?) ~ Ai{^i — fi). 
The complementary information on the non-local corre- 
lation functions of the bulk system is instead completely 
missing, and indeed the trapped system cannot faithfully 
reproduce the correlation properties of the bulk system 
in general (as testified by the poor performance of the 
LDA at the level of correlations [36|). Nonetheless, the 
trapped system can faithfully reproduce the correlations 
around the trap center over a length scale corresponding 
to the extent of the quasi-homogeneous central region. 
This aspect will be discussed in details in Section |VI[ 
where it will be shown in addition that the evolution 
of global correlation properties of the system (as probed 



e.g. by time-of- flight measurements) is dominated by the 
evolution of correlations in the trap center, so that cor- 
relations in the trap center can effectively be accessed in 
the experiments. 

The quasi-homogeneity condition Eq. ([t]) allows to 
identify nc as a close approximation to the ground-state 
density of the bulk Hamiltonian Eq. Q (with Vt — 0) 
at a chemical potential corresponding to the background 
chemical potential ^{Vt, N, J,U,V2) of the trapped sys- 
tem. As further elaborated in Section [Vj nc is experi- 
mentally accessible with conventional methods. In order 
to convert the information on nc into information on the 
bulk phase diagram of the Hamiltonian Hq of Eq. ^ , we 
need to know at which chemical potential the experiment 
is operating once the experimental parameters Vt, N, J, 
U and V2 are set. This is the goal of next Section. 



III. TRAP-SQUEEZING CONTROL OF THE 
CHEMICAL POTENTIAL 

In this section we investigate the dependence of the 
background chemical potential /Lt(Vt, TV, J, C/, V2), stabi- 
lizing a ground state with N particles, for the Bose- 
Hubbard Hamiltonian with parameters J, J7 in a trap- 
ping potential of strength Vt and, generally, in a superlat- 
tice potential of strength V2. We review the conventional 
Thomas-Fermi approximation and its generic prediction 
for the functional form of /i in the general case of a d- 
dimensional hypercubic lattice. Remarkably, a numerical 
investigation based on numerically exact quantum Monte 
Carlo shows that the Thomas-Fermi prediction, modified 
by quantum fluctuations, is rather accurate at least in the 
regime of either weak enough interaction or high enough 
density. As a result, the chemical potential turns out to 
be simply controlled in the experiments either via the 
control on the particle population N or the trap strength 
Vt, or on both. 



A. Atomic limit and Thomas-Fermi approximation 

A simple theoretical approach giving the relationship 
between the chemical potential and the other Hamil- 
tonian parameters for the Bose-Hubbard model is the 
atomic limit ( AL) , which consists in discarding the quan- 
tum kinetic term in the Hamiltonian and in solving for 
the diagonal part. Such an approximation is reasonable 
in the strongly interacting and strongly trapped limit 
U, VtR^ ^ J, where R is the radius of the atomic cloud. 
Minimizing the potential energy part with respect to the 
density 

■Hpot = ^ X! ni{ni-l)+V2 ^ gin^+Vt ^ m-fi ^ 

i i i i 

(8) 
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one finds 



li + U/2- V2 - Vt {n - r^f 
U 



(9) 



Imposing the condition N = Ui and passing from the 
lattice to the continuum formulation, one readily obtains 
the AL chemical potential for a d-dimensional system 



U 



'd + 2 r(d/2 + l) 



2 7r'i/2 

'0.82548... ([/ 7V)2/3 Vl'^ {d=l) 

= <( 0.79788... ([/ iV Ft)^^^ (^^ = 2) (10) 

0.81346... ([/ iV)2/5 Vt'^ {d^i) 

Notice that the superlattice term does not enter in this 
formula because it has been chosen so as to take symmet- 
ric values around zero, and therefore it vanishes upon 
spatial integration. 



On the other hand, a similar formula to Eq. ( 13 1 can 



be obtained in the weakly interacting case via a stan- 
dard Gross-Pitaevskii (GP) approach plus a Thomas- 
Fermi (TF) approximation |35]. Taking the mean-field 
approximation ~ ^'i, a| ~ on the normal ordered 
Hamiltonian Eq. ([I]) (where 4"^ is the condensate wave- 
function) one obtains the lattice GP energy functional 



Egp = -J^(4'**j -fc.c.) + 

iij) 



i i 

(11) 

where Vi = ¥2 gi + Vt {n - tq)^ - /x. 

The standard TF approximation consists in neglect- 
ing completely the kinetic term in the Gross-Pitaevskii 
functional. This is justified in the continuum case be- 
cause the kinetic energy is suppressed when the wave- 
function is slowly varying in space; on the other hand, 
in the lattice slowly varying wavefunction is such 

that vE'*^'j w l^'jpi namely the kinetic term does not 
cancel, but it effectively adds up to the chemical poten- 
tial term, — )■ /i — 2dJ. Minimizing the GP functional 
leads then to the lattice TF equation for the density: 



fi + 2dJ- V2 - Vt in - rpf 
U 



(12) 



Integrating over space in the continuum limit leads to the 
result ^TF — —2dJ + ^al + U /2. 

The central prediction is the linear dependence of /i on 
the combination (A^2y^d-)i/(2+d)^ with a slope dependent 
on both predictions will be verified by a nu- 

merically exact calculation in a large parameter regime. 
What is completely missing in the above simple approach 
are: 1) lattice commensuration effects, relevant in the 
case of the appearance of band/Mott insulator states 
(given that the analytical expression Eq. ( 10 ) is obtained 
in the continuum limit); 2) a proper treatment of the 
quantum kinetic J-term in the Hamiltonian. The lattice 
effects can be readily restored by numerically performing 



the sum N = J^i '^i instead of analytically integrating the 
continuum generalization of Eq. [9j at the expense of los- 
ing the closed-form prediction of Eq. ([To|). On the other 



hand, a full account of the quantum corrections requires 
a more extensive numerical treatment. This will be pro- 
vided in the following, where we will see that such effects 
do not alter too drastically the TF/AL predictions, and 
they can be captured already at the mean-field level. 



B. Quantum Monte Carlo results 

Here we present the results of quantum Monte Carlo 
simulations of the d-dimensional Bose-Hubbard model 
on L'' hypercubic lattices, based on the Stochastic Se- 
ries Expansion method with directed-loop updates |39|. 
We have performed simulations at low temperatures T ^ 
J/L in order to remove significant thermal effects. We 
perform the simulation in the grand-canonical ensemble, 
and we fine-tune the chemical potential /i which sta- 
bilizes a given particle number N for the Hamiltonian 
Eq. ([1]); this allows us to numerically sample the func- 
tion ii{Vt, N, J, U, V2). In the following we take J as the 
energy scale, and express all other quantities in units of 
J. We have performed the simulations for two values 
of the U/J ratio for each dimensionality d and both for 
zero superlattice potential (V2 = 0) and for an intense su- 
perlattice potential (V2 = U). We have considered both 
incommensurate superlattices (a = 830/1076 as in the 
experiment of Ref. I22p and commensurate ones (a = 3/4 
and 1/2). Results for one spatial dimension have already 
been partially reported in Ref. [501 



1. d=l without superlattice 

We begin our discussion with the case of absence of 
a superlattice (V2 — 0). Figs. [2] and Figs. [3] show re- 
sults for the two complementary regimes of U / J = 5 
and U/J = 20. Judging from the phase diagram of the 
Id Bose-Hubbard model 40 , for U / J = h the kinetic 
and the potential part of the bulk Hubbard Hamiltonian 
are in strong competition for most filling values n < 
so that a variation of the chemical potential brings the 
system through an alternation of correlated superfiuid 
phases and Mott insulating phases separated by quantum 
critical points. Despite the important quantum effects 
taking place in the system, it is remarkable to observe 
in Fig. [2] that all /i values obtained by various combina- 
tions of Vt and N collapse onto the same universal curve, 
which is a homogeneous function of N'^Vt/J, in agree- 
ment with the AL/TF prediction Eq. ([lO|. Refs. [4T1I42] 
have introduced the the so-called "characteristic den- 
sity" , p — N{Vt/ jy^"^, as the relevant parameter to char- 
acterize the density of a trapped system in d-dimensions. 
Our data suggests that the chemical potential is a ho- 
mogeneous function of the characteristic density. This 
result will be further confirmed in the case d = 2,3; see 
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also Sec. HID for a general discussion. 

Even more remarkably, for large enough N and/or 
Vt/J values, the universal curve obeys a modified TF 
(mTF) scaling of the kind 

tlmTF{VuN,J,U,V2)/J = 

Cd{U/J)^' +jd(U/J N)^^ iVt/J)^^ . (13) 



where Cd — Cd{U / J,V2/ J). Here 7^ is a d-dependent 
slope: numerical evidence shows that 7^ does not seem 
to depend on the other Hamiltonian parameters (see 
Ref.l5Dl and the discussion below). A fit to the d=l data 
for U/J = b gives us 7^=1 = 0.817(2), which compares 
surprisingly well with the AL/TF prediction of Eq. (10 1. 
Moreover Cd is a d-dependent offset term which is also 



found to depend on [//J and on V2 / J, and which contains 
the most relevant quantum corrections to the TF/AL re- 
sult. It is important to stress that in general Cd [U / J) 
is not the correct value of fi/J for Vt/J,N — >• 0. In fact 
in this limit the system becomes extremely dilute and 
the dependence on U/J should drop out. Indeed, in the 
limit of low N and Vt, the /i curve crosses over to the 
correct dilute limit (see inset of Fig. [2| which gives the 
well known result /i(iV = 0, = 0) = —2dJ as obtained 
from the solution of the tight-binding model. As seen in 
Sec. |III A this result coincides with the N = Q,Vt =0 



value predicted by the lattice TF theory (Eqs. ( 10 ) and 
( 12 )) - as it should be expected, given that, in the dilute 
limit, the bose gas becomes ideal. This also means that 
the quantum offset term (—2dJ) of the chemical poten- 
tial from lattice TF theory does not correspond at all to 
the quantum correction appearing in the mTF scaling, 
Eq. (13). Hence the form of Eq. (13 1 is far from being 
trivial. 




FIG. 2: Chemical potential for the Id Bose-Hubbard model 
with U = 5 J. The universal curve fi = ^{N'^Vt / J) has been 
obtained via the collapse of many data sets, with A'' = 20, 
...,120 and Vt = 0.004, 0.140 J. The solid curve is a linear 
fit to the Ansatz Eq. (13 1. Inset: zoom on the crossover from 



the mTF behavior to the low-A'^ and low-Vt behavior. 



The results for U = 20J are instead more elaborate. 
In this case the inspection of the Id Bose-Hubbard phase 
diagram reveals that the system with a filling n < 3 is 



a Mott insulator for most of the chemical potential val- 
ues. Hence the lattice stabilizes commensurate insulating 
regions in the trap which are completely missed in the 
continuum limit leading to the AL/TF formula Eq. (10 1. 
Nonetheless it is remarkable to observe that all the data 
obtained for different N and Vt values collapse onto the 
same universal curve which is again a homogeneous func- 
tion of N'^Vt, consistently with the AL/TF prediction. 
Yet this curve shows a succession of kinks separating var- 
ious regimes of filling and trapping: an inspection in the 
microscopic structure of the states corresponding to the 
various regions shows that each kink marks the appear- 
ance of particles in a new "shell" of the trapped cake 
structure of the density profile (see Fig. [3]ja)). When 
the system has filling (n^) < 1 throughout the trap, the 
large U /J value suppresses multiple occupation and the 
system behaves effectively as a hardcore-boson system. 
This is shown by the excellent agreement between the 
\ow-N /\ow-Vt data for the softcore boson system and the 
same data for a hardcore-boson system which can be 
obtained exactly by standard Jordan- Wigner diagonal- 
ization |43j . Increasing the filling and/or the trapping, 
the condition of single occupancy is eventually violated 
and the second shell of the cake starts to be filled. In 
the AL, the condition for the presence of a site with oc- 
cupation n = 2 in the center of the trap in a Id sys- 
tem is Vt (iV/2)^ > U which gives the critical value 
[N'^Vt/ J)n=2 — 4 U/J. Assuming that the population 
of further shells appears upon depletion of the shell with 
single occupation, we find that the n = 3 shell starts 
being occupied when 



Vt {Njlf = Vt {N2/2Y + U = 2U 



(14) 



which means that the n 
increasing Vt and the n 



2 shell stops growing upon 

3 shell starts building up. 
Here Ni indicates the number of particles in the i-th shell. 
With the condition Ni 



N2= N Eq. ( 14 1 gives 



(iWt/J)„=3 = 4(1 + V2f U/J 



(15) 



A similar reasoning leads to the critical value for the 
population of the n = 4 shell 

(iVVt/J)„=4 - 2(1 + V2fil + Vsf U/J (16) 

and so on. Hence for U = 20 J we obtain {N'^Vt/ J)n = 
80, 466.27, and 1740.16 for n = 2,3, and 4; the corre- 
sponding numerical estimates from the QMC simulations 
are {N'^Vt/J)n ~ 66, 400 and 1200 which, not surpris- 
ingly, is lower than the AL prediction due to quantum 
effects allowing the particles to tunnel from the wings of 
the cloud into the center at a lower trapping strength/ 
particle number. 

It is evident from Fig. [3] that lattice commensuration 
effects become weaker and weaker on the /i curve for in- 
creasing filling. In particular, when plotting the curve 
as a function of the natural TF parameter (N^Vt/JY^^, 
as done in Fig. [3f^c), we observe that the mTF scaling of 



7 



Eq. 13 sets in approximately once the n = 3 shell starts 



to be filled. A fit of the large- A^/large-Vt data with Eq. 13 
delivers a value of 7^=1 consistent with the value obtained 
from the U j J — 5 data, confirming that only depends 
on dimensionality. 
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enough filling follows the mTF behavior Eq. 13 and a fit 
to that equation gives again a 7^ consistent with what 
found in the previous section in absence of the superlat- 
tice. This confirms the TF/AL prediction that the pres- 
ence of a superlattice is not influencing the overall slope 
of the /i curve. Similar curves for the incommensurate 
superlattice case a = 1076/830 have already appeared in 
Ref.iSOi 




5 10 



5 10 



FIG. 3: Chemical potential for the Id Bose-Hubbard model 
with t/ = 20 J. a) Main panel: Universal curve = 
li[N'^Vt/J). The curve is obtained via the collapse of many 
data sets (same parameters as in Fig. [If. The differently col- 
ored regions correspond to the number of shells populated in 
the "wedding cake" structure, from 1 up to 4. Upper pan- 
els: four representative density profiles showing four different 
shell occupations (iV = 100 and Vt = 0.004, 0.014, 0.070 and 
0.135 J), b) Zoom on the low- TV/low- Vi region (first shell oc- 
cupied only), showing agreement with the corresponding data 
for a hardcore boson system, c) Universal curve presented as 
a function of N^^^{Vt/ Jf'^: the solid curve represents a lin- 



ear fit to the Ansatz Eq. ( 13 1 



2. d=l with superlattice 

Fig. |4] shows the fx curve for the Id Bose-Hubbard 
model in a strong superlattice with a — 3/4, (j) — Q and 
V2 = C/, and for three values of U/J = 10, 20 and 30. 
Again the chemical potential is clearly a homogeneous 
function of N'^Vt/ J , as shown by the collapse of all the /i 
values onto the same universal curve dependent on U / J 
and V2/J. For large U/J, V2/J and/or for smaU N'^Vt/J 
we observe kinks in the /i curve characteristic of lattice 
commensuration effects: in the case of a commensurate 
superlattice as the one considered here these effects are 
clearly related to the fractional filling shells [30] appear- 
ing in the AL of the system at fillings (2n -|- l)/4 with 
n = 0, 1, 2, .... Yet the overall trend of the curve for high 



FIG. 4: Chemical potential for the trapped Id Bose-Hubbard 
model with a commensurate superlattice potential (V2 = 
U ,a — 3/4). The right panel shows that all curves have a 
universal scaling of the the slope with the ratio U /J. Dashed 
lines are fits to the Ansatz Eq. (13 1. 



It is important to stress that lattice commensuration 
effects in the presence of a superlattice are strongly re- 
lated to the possibility of maintaining the spatial phase (j) 
of the superlattice fixed in the experiments. If this phase 
is allowed to fiuctuate, (as it generally happens from shot 
to shot of the same experiment [H] , unless phase locking 
is explicitly enforced [H]), the chemical potential in the 
experiment is going to change accordingly, assuming that 
all other experimental parameters remain unchanged. In 
this case it is then more convenient to introduce a phase- 
averaged chemical potential (fJ-)^ - namely averaged over 
random fluctuations of the phase </>. Fig. [5| shows (/i)^ for 
two superlattices: for the commensurate one discussed 
above, and for an incommensurate one which has one 
additional color component, namely 

gi{a, a' , 0, 0') = cos^(27rQ; i + (p) + cos^(27ra' i + (f)') — 1 

(17) 

with a = 1076/830 and a' = 1473/830. This particular 
superlattice will be further discussed in Section [VlT| The 
chemical potential values are typically averaged over a 
sample of ~ 100 — 200 phase values. Comparing the av- 
erage results with typical ones we notice that fluctuations 
are small around the average, so that, for each different 
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U/J 


V2 


a 


7d=l 




5 





- 


0.818(1) 


-1.23(1) 


10 




- 


0.818(2) 


-1.21(2) 


20 




- 


0.815(5) 


-1.32(4) 


10 


u 


3/4 


0.817(2) 


-1.23(2) 


20 




3/4 


0.820(2) 


-1.40(1) 


30 




3/4 


0.822(3) 


-1.58(1) 


10 


u 


830/1076 


0.816(3) 


-1.22(3) 


20 




830/1076 


0.822(4) 


-1.42(4) 


30 




830/1076 


0.819(5) 


-1.54(5) 



TABLE I: Results of the fit of the chemical potential data 
to the scaling Ansatz, Eq. ( |13| l, for the Id trapped Bose- 
Hubbard model without any superlattice, and with a com- 
mensurate / incommensurate superlattice. 
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realization of the superlattice, the average chemical po- 
tential (n)^ stabihzes in the system a number of particles 
N which is close to the desired one. 



FIG. 6; Chemical potential for the trapped 2d Bose-Hubbard 
model without (left panel) and with (right panel) a super- 
lattice potential (a = 830/1076 in all spatial dimensions). 
Dashed lines are fits to the Ansatz Eq. (13 1. 



60 p 
50- 

40- 



A 30 
^ 20 
10 

0, 




o commensurate superlattice 

□ incommensurate 2-color superlattice 



6 7 8 9 
1/3 

N- (V /J) 



10 



FIG. 5: Phase-averaged chemical potential for the trapped 
Id Bose-Hubbard model with a commensurate superlattice 
potential and with a 2-color incommensurate potential (see 
text). For both cases V2 = U ^ 20 J and N = 100. 



As a summary of the study of one dimensional sys- 
tems, Table |I] shows the result of the fit of the various 
cases considered (no superlattice, commensurate and in- 
commensurate superlattice) to the Ansatz of Eq. (13 1. 



We observe that the coefficients jd^i are essentially all 
consistent within error bars, as the TF/AL theory would 
simply predict; more quantitatively, they are all close to 
the TF/AL prediction 7^=1 = 0.82548.., which is a re- 
markable fact given the simplicity of that theory. More- 
over the constant terms Cd=i{U/ J,V2/ J) appear to be 
only weakly dependent on the Hamiltonian parameters; 
in particular, within error bars they appear not to de- 
pend on the parameter a of the superlattice but only on 
the superlattice strength. 



3. d=2, 3 

Fig. |6] shows the chemical potential for the d — 2 Bose- 
Hubbard model as numerically determined via quan- 
tum Monte Carlo. We have considered the two val- 
ues U/J= 15 and 40 for the interaction: according to 
the phase diagram of the homogeneous 2d Bose-Hubbard 
model [45], for the first value the system is always in a 
superfluid state at all fillings, while for the second value 
the system experiences Mott insulating phases at fillings 
n = 1,2, .. upon changing the chemical potential. From 
Fig. [6] we see that, no matter the phase the system is 
in, data sets for different particle numbers and varying 
trapping potentials all collapse onto the same curve when 
plotted as a function of {NVtY^"^, agreement with the 
TF/AL prediction, both in absence and in presence of 
a superlattice potential. For sufficiently high N and/or 
Vt the dependence of /i on {NVtY^'^ linear and can 
be well fitted to the Ansatz Eq. ( 13 1 giving the results 
summarized in Table [Hi It is remarkable that the nu- 
merical results for the slope 7^=2 are very close to be- 
ing all numerically consistent with each other, and with 
the TF/AL prediction 7^=2 = 0.79788. On the oppo- 
site end, for small N and/or low Vt, the local filling can 
drop to values {rii) < 1, which, in absence of a super- 
lattice, amounts to the onset of a n = 1 Mott plateau 
for U/J = 40: this strong lattice commensurability ef- 
fect, completely neglected in the TF/AL approach in 
continuum space, is responsible for the deviation of the 
U/J = 40 results from the linear behavior. As already 
observed for the d ~ 1 case [30] the fj. curve exhibits a 
crossover to a low-density regime which loses the depen- 
dence on the U/J value, as shown by the merging of the 
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7d=2 
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0.795(3) 


-2.1(1) 


40 
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0.799(3) 


-3.53(8) 



TABLE II: Results of t he fi t of the chemical potential data to 
the scaling Ansatz Eq. ( 13 1 for the 2d trapped Bose-Hubbard 



model without any superlattice, and with an incommensurate 
superlattice. 
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-4.55(11) 
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-4.4(2) 


50 
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0.81(5) 
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TABLE III: Same as in Table |TT] bur for the 3d trapped Bose- 
Hubbard model. 



U/J — 40 curve with the U/J — 15 one. In general one 
does expect all fi curves at arbitrary U/J values to merge 
onto the same universal curve for a diluted system when 
N,Vt — )• 0. Yet a large enough U/J suppresses double 
occupancy in a system with filling smaller than unity, 
so that the the actual value of U /J becomes irrelevant 
and the system enters the so-called Tonks or hardcore 
regime: this is the effect responsible for the merging of 
the fj. curves at U/J = 40 and C//J = 15 in Fig. |6j well 
before the regime of extreme dilution is attained (the 
two curves come close to each other for (NVtY^^ values 
which correspond to a filling n ~ 1 in the trap center). 
On the contrary, for the NVt values we considered, the 
hardcore regime is not observed in presence of a strong 
superlattice V2 = U, which competes with the repulsion 
and maintains doubly occupied sites down to the lowest 
trap fillings we have explored. 

To conclude. Fig. [7] shows analogous results for the 
d — 3 case. We again choose two values of repulsion: 
U/J = 20, for which the homogeneous system is super- 
fluid at all fillings [15]; and U/J = 50, for which the 
system has Mott insulating regions at fillings n — 1,2... 
. In agreement with the TF/AL prediction, the chemi- 
cal potential is a homogeneous function of the product 
iV2/5y^3/5^ and it becomes linear for large enough filling 
in the trap center (n > 1); the results of fits to Eq. (13) 



are summarized in Table |III[ and show numerical consis- 
tency between the estimated slopes 7^=3 and the prox- 
imity to the TF/AL prediction 7^=3 = 0.81346, both in 
absence and in presence of a superlattice |47) . In the op- 
posite limit of low filling n < 1, and in absence of the su- 
perlattice, we observe again that the /i curves at different 
U/J values tend to merge together, which corresponds to 
the onset of the hardcore regime with suppression of dou- 
ble occupancy. 
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FIG. 7: Chemical potential for the trapped 3d Bose-Hubbard 
model without (left panel) and with (right panel) a super- 
lattice potential (a = 830/1076 in all spatial dimensions). 
Dashed lines are fits to the Ansatz Eq. (13 1. 



Quantum Monte Carlo vs. Gutzwiller 
mean-field results 



In the previous section we have seen that the chemi- 
cal potential of strongly correlated bosonic systems can 
be described via a modified Thomas-Fermi scaling for 
sufficiently high densities or sufficiently strong trapping 
potentials, at which lattice commensuration effects be- 
come very weak. On the other hand, lattice effects 
can be easily taken into account at the classical level 
via a numerical calculation in the atomic limit, neglect- 
ing intersite hopping. Yet such a drastic approximation 
is going to be unreliable when the hopping energy be- 
comes comparable with the repulsion one, typically when 
2dnJ ^ Un?. Nonetheless, quantum effects can be re- 
stored at an approximate level within Gutzwiller mean- 
field theory, which keeps the cost of numerical calcula- 
tions at a minimum, while producing predictions which 
turn out to be in a suprisingly good agreement with those 
of quantum Monte Carlo. 

Gutzwiller mean-field theory assumes a factorized 
Ansatz for the wavefunction of the many-body system 




(18) 



where \n)i is a Fock state with n particles at site i, and 
f^max is a suitable truncation imposed on the local Hilbert 
space for numerical purposes. Assuming all fn^ coeffi- 
cients to be real, the expectation value of the Hamilto- 
nian, Eq. ([T]), on the Gutzwiller Ansatz (GA) wavefunc- 
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FIG. 8: Chemical potential for the trapped Id Bose Hubbard 
model with U = lOJ, as obtained via quantum Monte Carlo, 
Gutzwiller Ansatz and atomic limit calculations. /XmTF is the 
modified Thomas-Fermi prediction of Eq. ( |13[ ) with fitting 
parameters from Table |l] Upper panel: V2 = 0. Lower panel: 
V2 = U,a^ 830/1076. 



tion takes the form 



(19) 



to the iVriinax coefScients. Instead of proceeding with a 
full minimization, a usual procedure is to minimize with 
respect to the local variables /i,*-* at a site i while hold- 
ing fixed the variables at all the other sites. We perform 
random sweeps over the lattice sites, touching each site 
once per sweep on average, until convergence is reached 
- typically less than a hundred sweeps are necessary for 
convergence. For the Hubbard model under investiga- 
tion, this leads to a very efficient location of the absolute 
energy minimum, which can be verified by repeating the 
minimization starting from a different initial condition. 
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where 



and 



— n(n — 1) + Viii 



I n 



V2 gi + Vt{ri - r^f - ^ 



(20) 



(21) 



Here runs on the z = 2d nearest neighbors in a hy- 
percubic lattice. 

The ground state corresponds then to the minimiza- 
tion of the Hamiltonian expectation value with respect 



FIG. 9: Chemical potential for the trapped 2d and 3d Bose 
Hubbard model. /^mTF is the modified Thomas-Fermi predic- 
tion of Eq. ( 13 1 with fitting parameters from Tables |ll] and 
|III[ Upper panel: d — 2, U — 15J, V2 ~ 0. Lower panel: 
d = 3,U = 20J, V2 = 0. 



Eq. ( 18 1 generally describes a state which does not 



have a well defined particle number. Given that the min- 
imization of the energy is done in an unconstrained way 
at each site, we need to adjust the chemical potential /i 
a posteriori in order to achieve a desired average particle 
number N — J2i „ Ifn^'l'^n, with an analogous procedure 
to that used with quantum Monte Carlo. Hence this 
procedure samples the function fi = ^{Vt,N, J, U, V2) of 
interest. We use different particle numbers N and trap- 
ping potentials, ranging in the same intervals as those 
explored in quantum Monte Carlo calculations. 
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Figs. [8] and |9] compare the results of the GA with those 
of quantum Monte Carlo calculations for the trapped 
Bose-Hubbard model in dimensions d—1,2 and 3. Data 
obtained in the atomic limit J = are also shown as a 
reference. All data are normalized to the mTF behavior, 
Eq. (13), attained in the large density / strong confine- 
ment limit. The results in the atomic limit are seen to 
generally overestimate the chemical potential - this is not 
surprising, given that there is no kinetic energy gain in 
adding particles to the system in the atomic limit, and 
hence all the energy gain to balance repulsion has to be 
provided by the chemical potential. Quite surprisingly, 
instead, we observe that the mean-field GA accounts very 
well for the quantum corrections to the atomic limit for 
all system dimensions, and that the chemical potential 
predicted via the GA is typically less than 4% off with 
respect to the quantum Monte Carlo prediction. Indeed 
the GA, while successfully describing the main features 
of the phase diagram of the Bose-Hubbard model, is not 
particularly accurate in determining the phase bound- 
aries; given its mean-field nature, it is not quantitatively 
trustworthy in low dimensions, and its predictions on the 
location of phase boundaries can be up to 100% off in 
d — 1. We find that the GA works equally well even 
in presence of an incommensurate superlattice potential, 
despite its inability to describe Anderson localization and 
Bose-glass physics (data for the case d = 1 are shown in 
Fig.|8]). 

A semi-quantitative justification of the success of the 
GA in predicting the chemical potential of the trapped 
Bose-Hubbard model can be provided based on LDA. 
LDA predicts that the total density of the system can 
be written as 



N 



(22) 



where Hi = ^ — Vt(ri — tq)^ is the local chemical poten- 
tial, and n(/u) is the density as a function of the chemical 
potential for the bulk system. The GA might provide a 
very poor prediction nQfi^(fi) for the n{fi) function of the 
bulk system, especially in d = 1, and this might suggest 
that the global chemical potential /x which stabilizes N 
particles in the system should be also poorly estimated. 
Yet what matters in the determination of TV is not the 
whole Ti(/i) curve, but only its integral. If the difference 
An(/i) = n(/i) — nGAifJ-) oscillates from positive to neg- 
ative and back, it will be averaged to zero upon integra- 
tion. This is most likely the case for the Hamiltonian pa- 
rameters we have considered. As observed e.g. in Fig. |8j 
the most serious problems arise for low N, which, in a 
trapped system, implies a chemical potential excursion 
(from the tails to the center) over which the An function 
has completed only a few oscillations (if any) . 



D. Discussion 

From the above results we can conclude that the chem- 
ical potential for the Bose-Hubbard model in a trap and 
in a superlattice potential has a simple scaling form as 
a function of the experimental parameters. On the one 
hand an accurate determination relies in general on nu- 
merics due to the strong interactions in the system; yet 
we find that the chemical potential appears as a homo- 



geneous function of Ar2/(2+d) 



d/{2+d) 



and in partic- 



ular a linear function thereof for sufficiently high fill- 
ing in the trap center (typically n > 1 in absence of 
a superlattice, and even lower in presence of a super- 
lattice), and/or weak enough interaction, verifying this 
way a quantum-modified version of the TF/AL predic- 
tion, Eq. (13). Ref. [35] has shown that, under the as- 



sumption that LDA holds, one can prove the relationship 
p — Id{p/ J',U/ J), where p is the characteristic density, 
already introduced above, and Id is an unknown func- 
tion. By inverting the previous relation, this amounts 
to say that p is a homogeneous function of p, in agree- 
ment with our findings. Therefore our results allow to 
explicitly reconstruct the relation between p and p. This 
correspondence of our results with the LDA prediction 
suggests that p is exactly a homogeneous function of p if 
(and possibly only if) LDA is exact. Indeed deviations 
from perfect homogeneity are seen in our numerical data, 
for instance in Fig. [8] and they can be attributed in part 
to numerical uncertainty; in part to the fact that in any 
finite system different p values correspond to the same 
N value (which is discretized), so that there is a natu- 
ral uncertainty on p; and in part to possible systematic 
deviations, enhanced in the presence of a rapidly oscil- 
lating superlattice. Yet we observe that the homogeneity 
property is an important simplifying assumption, and it 
is verified with an accuracy of a few percent in the worst 
case. 

Moreover, in the case of the Bose-Hubbard model with- 
out superlattices, Ref. [13] has numerically shown that 
density profiles (and profiles of other local observables) 
for trapped bosons having the same characteristic den- 
sity p, and the same Hamiltonian parameters U and J, 
appear to be invariant up to a rescaling of the space coor- 
dinates with the characteristic length rj = ^ J/Vt. This 
implies that, under rescaling of the coordinates with ry, 
the density at site i, (n;) = {h{ri/ri)), is a unique func- 
tion of p. At the same time, the local chemical potential. 
Pi, expressed in terms of rescaled variables, turns out as 
well to be a unique function of p 



= p{p-U,J) + J[{n~ra)/Tf 



(23) 



thanks to the result that p = p{p;U,J). As a conse- 
quence, the results of Ref. [49] suggest that {fi{ri/-q)) is 
a unique function of Pn/rj- In the limit of an infinitely 
shallow trap, Vt — and — >■ oo at fixed p, the de- 
pendence of {h(ri/ri)) on Pn/rj must reproduce the de- 
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pendence of the bulk density on the chemical potential, 
n{fi). Hence the unique function relating {h{ri/ri)) to 
ij-n/i-i must be identical with that of the bulk limit. As a 
consequence, the invariance of the density profiles (and, 
with the same reasoning, of any other local observable in 
the trap) under rescaling with 77 implies that the LDA 
is exact. It is trivial to show that the the exactness of 
the LDA implies the invariance of local quantities under 
rescaling with rj. Hence we can conclude that the scale 
invariance of trapped systems at fixed characteristic den- 
sity is exact if and only if LDA is exact. Yet, even when 
LDA is not exact, the system might display scale invari- 
ance with a better accuracy than that of LDA, as shown 
numerically in Ref. |49j in the case of density fluctuations 
and local compressibility. 

At the experimental level, the above results reveal 
the possibility of controlling the chemical potential in 
the trapped system by: 1) controlling the trapping fre- 
quency at fixed particle number, or viceversa 2) control- 
ling the particle number at fixed trap strength, or 3) 
controlling both simultaneously. Hence the present ex- 
perimental setups have direct access to the measurement 
of phase diagrams of strongly correlated bosons in the 
grand- canonical ensemble at variable chemical potential. 
Most noticeably, an accurate knowledge (within a few 
percent) of the zero-temperature chemical potential as a 
function of the Hubbard Hamiltonian parameters can be 
obtained with little numerical effort via the mean-field 
Gutzwiller Ansatz, namely without fully solving numer- 
ically the many-body problem. This aspect is to be con- 
trasted with the case of accurate thermometry of strongly 
correlated bosons, which at the moment can be reliably 
assessed only via exhaustive ab-initio simulations [T71ll8j . 



varying the trap strength Vt independently of the Hamil- 
tonian parameters U, J, V2, offering in this way the pos- 
sibility of simulating the Bose-Hubbard Hamiltonian in 
a variable trapping potential [5T1 [5^ . 




FIG. 10: Schematic view of the confinement control on the 
longitudinal direction in a set of one-dimensional tubes. 

In the remainder of the paper we hence focus on the 
proposal of extracting the phase diagram of the Bose- 
Hubbard model in the grand-canonical ensemble by con- 
trolling the chemical potential of the system via trap- 
squeezing and at fixed particle number. Two fundamen- 
tal issues are addressed in the next subsections concern- 
ing the experimental feasibility of this proposal: 1) we 
discuss the fundamental difficulty in achieving adiabatic 
trap-squeezing once the particles are loaded in the opti- 
cal lattice; 2) hence we propose a simple loading sequence 
which overcomes this difficulty, and which allows to con- 
trol the particle number in the central tubes/layers in the 
case of Id/ 2d confinement. 



IV. TRAP-SQUEEZING PROTOCOL 

The previous section has shown how the chemical po- 
tential of a trapped system is simply controlled by the 
number of particles N and by the trapping potential Vt- 
Indeed recent experiments have probed the compress- 
ibility of fermions loaded in optical lattices by varying 
the trapping potential [3T] or the atom number [SUl EI] , 
although the explicit relation between the experimental 
parameter and the chemical potential was not known for 
that system. In the experiments, the atom number is typ- 
ically subject to significant fluctuations due to shot noise 
of order \/N, as well as to other systematic effects. On 
the contrary the trapping potential Vt can be controlled 
on a much finer scale via the application of a dipolar 
trap along the spatial direction of interest, namely by 
shining d red-detuned running laser waves onto the sam- 
ple to control the confinement in d spatial dimensions 
(this is sketched in Fig. 10 for the case of variable one- 
dimensional confinement). The dipolar trap adds up to 
the confining (anti-confining) potential given by the over- 
all gaussian intensity profile of the red-detuned (blue- 
detuned) optical lattice, and it gives the possibility of 



A. The adiabaticity issue 

Optical lattice experiments have revealed the unique 
possibility of tuning the Hamiltonian parameters in real 
time, e.g. via changing the intensity of the standing wave 
which controls the U/J parameter |53j . A fundamental 
issue raised by such real-time control is the possible gen- 
eration of excitations in the system by a non-adiabatic 
change of the Hamiltonian parameters. This becomes 
particularly dramatic when the parameter change implies 
the crossing of a quantum critical point, at which the 
gap over the ground state vanishes, leading to efficient 
Landau-Zener tunneling j54j . Similarly to the change in 
the optical lattice strength, a change in the dipolar trap- 
ping strength can also lead to quantum phase transitions, 
driven this time by the chemical potential, and hence it 
is subject to adiabaticity issues. In fact adiabaticity is 
even more problematic in the case of trap squeezing: re- 
garding the center of the trap as the system, and the 
trap periphery as the particle reservoir (as mentioned in 
the Introduction), the process of transfer of particles be- 
tween system and reservoir under trap squeezing might 
be pathologically slow if the particles are strongly local- 
ized in one part or the other, due to strong interactions 
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in presence of an optical lattice or to (quasi-) disorder 
potentials. 
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FIG. 11: Density profile as a function of trap strength for 
the Bose-Hubbard model with A*' = 100 particles and with 
U = 20/ J, as resulting from an extended fermionization cal- 
culation (see text and Appendix [A|) . 
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FIG. 12: Lowest energy gap A as a function of trap strength 
for the Bose-Hubbard model. Parameters as in Fig. |11| 



We exemplify the above discussion in the case of lattice 
hardcore bosons, for which the exact many-body spec- 
trum can be calculated exactly via Jordan- Wigner diag- 
onalization 03], in the case of density n < 1 on each site. 
The fermionization approach can be extended to the case 
n > 1 |55j under the condition of having a well developed 
layer-cake structure of the density profile (see Appendix 
[A| for a detailed discussion). Fig. 11 shows the evolu- 



tion of the density profile upon trap squeezing for the 
trapped Bose-Hubbard model with U — 20J and with 
N — 100 particles, as obtained via extended fermioniza- 
tion. For the chosen range of trap strengths, we observe 



that the system goes from a Mott-insulating phase with 
n = 1 particles in the trap center to a superfluid phase 
with 1 < n < 2 and finally to a Mott-insulating phase 
with n — 2 particles. In particular, the "local transition" 
from Mott insulator to superfluid in the trap center at 
fixed particle number is exclusively due to the redistri- 
bution of particles from states in the wings of the trap 
to states in the center. The states in the trap periph- 
ery are localized by the joint effect of the strong repul- 
sion exerted by the particles in the trap center and of 
the confining potential. On the other hand the states 
in the trap center are localized by the confining poten- 
tial over a length scale which is well below the width 
of the atomic cloud. Hence the spatial overlap between 
such states is negligible, and they are separated by a 
high potential barrier - which is represented by the atoms 
forming the n = 1 Mott plateaus. This implies that the 
tunnel splitting between the states in the center and the 
states in the wings can be extremely small, and hence 
the ground state can become nearly degenerate with the 
first excited state. This is indeed revealed by the di- 
rect investigation of the lowest energy gap, which, within 
the extended fermionization approach, is estimated as 
the lowest particle-hole excitation energy for the spinless 



fermions, and it is shown in Fig. 12 We observe that such 
a gap goes to a value very close to zero 56] when the first 
particles move from the trap wings to the trap center, cor- 
responding to the occurrence of an insulator-to-superfluid 
transition in the trap core. The gap shows as many dips 
as particles transfered to the center, revealing the tightly 
avoided level crossings corresponding to the particle re- 
distribution - more precisely, each deep minimum of the 
A vs. Vt curve corresponds to the successive migration of 
two particles from the wings to the center, because such 
particles occupy states which are nearly degenerate sym- 
metric/antisymmetric superpositions of localized states 
on the two opposite wings of the trap) . Hence this suc- 
cession of tightly avoided level crossings makes adiabatic 
trap squeezing essentially impossible in presence of Mott- 
insulating regions, as squeezing times which are several 
orders of magnitude bigger than the typical tunneling 
time (~ms) would be required. 

As a second example we consider the case of A'' = 50 
hardcore bosons in an incommensurate superlattice po- 
tential Eq. ([3]) with a = 830/1076 (as experimentally 
realized in Ref. [22]) and V2 = 20J. Such a potential 
is known to lead to Anderson localization of all single- 
particle states for V2 > 4 J [371 ISH] ■ In this case standard 
fermionization provides the exact spectrum, and in par- 
ticular the lowest energy gap as a function of the trap 
strength, as shown in Fig. [T3j We observe that the gap 
can become extremely small, this time due to the quasi- 
degeneracy of Anderson localized states which are spa- 
tially well separated, so that the tunnel splitting between 
them is very small. 

Hence from these two examples we can conclude that 
either strong repulsion induced by the optical lattice, or a 
(pseudo-) disorder potential, give rise to localized single- 
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FIG. 14; Sketch of the proposed loading and squeezing pro- 
FIG. 13: Lowest energy gap A as a function of trap strength ^^^^^ ^^^^^^ spatially anisotropic optical lattice. The 
for iV = 50 hardcore bosons in an incommensurate superlat- description of the various stages is presented in SeclrVBl 
tice (q = 830/1076) of strength V2 = lOJ. c I r 



particle states, which become quasi-degenerate upon 
changing the trap frequency. This leads to a very ef- 
ficient Landau-Zener tunneling which essentially makes 
adiabatic trap squeezing impossible. Although the re- 
sults shown here are obtained for a one-dimensional sys- 
tem, their underlying mechanism is quite general and ap- 
plies to higher dimensions as well. Hence we can gener- 
ally conclude that adiabatic trap-squeezing for a system 
of particles loaded in a strong optical (super)lattice is 
virtually impossible. Luckily the statement of the prob- 
lem already contains in itself the solution: trap squeez- 
ing has to be performed before loading the particles in 
the optical (super)lattice, namely when the particles still 
enjoy their full mobility, so that they can adiabatically 
follow the variation in the trapping potential. In the case 
of spatially anisotropic optical lattices, discussed in the 
next subsection, this consideration applies to the load- 
ing of the particles in the optical lattice along the spatial 
dimensions in which the particles are more weakly con- 
fined. 



B. Trap squeezing protocol for low-dimensional 
systems 

In the case of spatially anisotropic optical lattices, al- 
lowing to simulate the Id or 2c? Bose-Hubbard model 
pn the control of the chemical potential can be 
achieved by trap squeezing along the spatial dimensions 
of interest, namely along the tubes in lc?-anisotropic op- 
tical lattices and along the layers in 2d-anisotropic ones. 
Yet, if the strength of the trapping potential is changed 
before confining the atoms in tubes or in layers, the final 
atomic population present in the tubes/layers will be sig- 
nificantly affected, so that one violates the condition of 
working at fixed particle number in the low-dimensional 
elements (tube or layer) of the system, and in partic- 



ular in the central ones. The problem of guaranteeing 
at the same time the control on the particle number in 
the tubes/layer and adiabatic trap squeezing can be eas- 
ily solved by considering the following loading protocol, 
sketched in Fig. [14} 

a) Particles are initially trapped in a dipolar trap or a 
magneto-optical trap (MOT) at fixed strength; 

b) A deep one-dimensional/two-dimensional optical 
lattice is then ramped up, defining layers/tubes in which 
the atoms are still moving in continuum space. For a 
given initial atom number N the overall trapping po- 
tential (dipolar trap or MOT plus the optical lattice 
beam profile) uniquely defines the populations in each 
tube/layer. These populations remain essentially fixed 
for the rest of the experiment, thanks to the deep optical 
lattice. Given that the atoms are still in the weakly in- 
teracting regime, we can determine these populations via 
the Thomas-Fermi approximation - the results are given 
by Eqs. (|B6]) and ([^ in Appendix [b]); 

c) Varying the dipolar trapping along the tube/layer 
dimensions gives rise to trap squeezing, which can be 
easily kept adiabatic for weakly interacting particles [30] ; 

d) At this stage the optical lattice can be adiabati- 
cally ramped up along the tubes/layers, finally realizing 
a ld/2d Bose-Hubbard Hamiltonian at a given particle 
number and in a given longitudinal trapping potential. 



V. IMAGING TECHNIQUES 

The two previous sections have discussed how to 
achieve full control on the chemical potential of the Bose- 
Hubbard model realized in optical lattice experiments via 
the control on the trapping potential at fixed particle 
number. As discussed in Section [H] the overall chemical 
potential of the trapped system corresponds in partic- 
ular to the local chemical potential at the trap center: 
hence, getting access to the average central density in 
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the trap as a function of the chemical potential gives the 
possibility of extracting experimentally the density curve 
of the bulk Bose-Hubbard model in the grand-canonical 
ensemble and hence its compressibility. 
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FIG. 15: Sketch of the proposed imaging protocol for alkali 
atoms. A repumping {UI2) and a cooling (wi) beam are tightly 
focused on the center of the trap, and they are crossed so as 
to overlap in a region of width ~ W (minimal beam waist) in 
all spatial directions. This setup realizes a spatially-resolved 
optical cycling transition. 



As it will be further seen in the next Sections, measur- 
ing an average central density which mimics closely the 
one of the bulk system requires imaging the atomic cloud 
over a typical length scale of 10-20 lattice sites in each 
spatial direction for typical atom numbers and trapping 
frequencies used in current experiments, corresponding 
to a spatial resolution of ~ 5 — 10 /xm. To achieve this 
resolution, a possible technique is given by spatially re- 
solved fluorescence imaging for atoms (see Fig. 15 1. In 
the case of alkali atoms initially prepared in an F 



Ihy- 

perfine state of the level, a repumping beam, reso- 

nant with the transition to an F' = 2 state of the n^P3/2 
level, optically pumps the atoms into an F = 2 state [61] ; 
a second beam resonant with the F = 2 ^ F' — 3 tran- 
sition (cooling beam) is then used to image the atoms. 
If the two beams are tightly focused to a beam waist of 
^5 — 10 fim and are perpendicular to each other such 
that they meet at the focal point, they selectively image 
atoms in a region of space of the desired size. Typical 
optical lattice experiments have n > 1 particles per lat- 
tice site in the trap center, so that imaging a central 
region with lO'^ — 20"^ lattice sites involves imaging at 
least as many atoms, which is possible over a time < 100 
ms (see e.g. Ref. Before imaging, the optical lat- 

tice could be rapidly increased to the maximum height 
in order to freeze the atomic cloud profile during the suc- 
cessive imaging time. 

The proposed imaging technique has the advantage of 
being generally compatible with common atom trapping 
setups based on magneto-optical or dipolar traps, and 
to require only an intermediate imaging resolution. On 
the other hand, ultra-high resolution optics has given ac- 
cess to few-site/few- atom or even single-site/single-atom 
imaging, in a variety of very recent experiments |25H29j . 
This level of resolution far exceeds the one required by 



our present proposal. An alternative technique for the in- 
direct extraction of the local density in the trap consists 
in the high resolution in-situ imaging of the atomic cloud 
density integrated along the line of sight, from which 
the full three-dimensional atomic distribution is recon- 
structed via inverse Abel transformation [63l Elj . Given 
the diversity of techniques mentioned above, we are confi- 
dent that the extraction of the local density properties in 
selected regions of the trap will become an experimental 
routine in the near future. 

After the in situ imaging stage, turning off all trapping 
potentials and imaging the expanded cloud gives access 
to the total number of atoms. This piece of information 
is fundamental to post-select the measurements with a 
given total atom number A^, which corresponds to the 
desired value of the chemical potential to be realized in 
the experiment. Strictly speaking, once the calibration 
curve relating /i to A is known, the final measurement 
of the particle number has only the role of assigning the 
measurement of the central density to the proper place 
in the grand-canonical phase diagram of the bulk system. 
Hence one can regard the fluctuations in the total parti- 
cle number TV as a source of random sampling of the fi 
axis in that phase diagram. In this perspective, all mea- 
surements give useful information, provided that particle 
number fluctuations are not bringing the chemical po- 
tential fi too far from the region of interest in the phase 
diagram. 



VI. APPLICATIONS: PHASE DIAGRAM OF 
THE Id BOSE-HUBBARD MODEL 

As a first application of trap squeezing, we investigate 
the ground-state phase diagram of the one-dimensional 
Bose-Hubbard model. We consider the situation of A = 
100 bosons trapped in the central tube of a strongly 
anisotropic optical lattice [3T], and addressed with spa- 
tially resolved imaging described in the previous section. 
We present quantum Monte Carlo results for the central 
properties in the trap, as well as for the global time-of- 
fiight properties. The chemical potential ^ is adjusted so 
that the ground state of the grand-Hamiltonian, Eq. |4j 
has the desired number of bosons. 



A. Central compressibility 

Investigating the phase diagram of the Hubbard model 
via trap squeezing implies that the superfluid-insulator 
transitions that one can probe are accessed via a vari- 
ation of the chemical potential, namely they are of the 
commensurate-incommensurate kind in a system with- 
out disorder. Such transitions have the advantage of be- 
ing very sharp, as they are characterized roughly speak- 
ing by the appearance (or disappearance) of a superfluid 
component in the system, induced by doping particles or 
holes into a Mott insulating state at integer filling. The 
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-^/^ K^U, 10 sites -^/^ K^U, 20 sites 

FIG. 16: Central compressibility = dn^/dfj, (in units of oi N = 100 bosons under variable trapping frequency. The 

central density is obtained by imaging over 10 sites (left panel) and 20 sites (right panel). The white curves report the 
boundary lines between superfluid and Mott-insulating phases as determined in Ref. |40| (for the n = 1 and n — 2 lobe) via 
density-matrix renormalization group, and independently in this work (for the n — 3 lobe) via quantum Monte Carlo. 



superfluid-insulator transition can be completely charac- 
terized by the compressibility k — dn/dfi, which, simi- 
larly to the superfluid density, is identically zero in the 
Mott-insulating phase and finite in the superfluid phase, 
for all spatial dimensions. In particular the compress- 
ibility has quite a sharp variation at the transition in 
low-dimensional systems. Indeed its critical scaling with 
system size L obeys the law L*^"^ [5] where z = 2 is the 
dynamical critical exponent; hence, for d = 1, it jumps 
from a value diverging linearly with system size to zero at 
the transition point (while in d = 2 the jump takes place 
from a value diverging logarithmically with the system 
size). This very sharp feature makes it possible to de- 
tect the transition unambiguously even on a small system 
size, which is effectively the case here. In the following 
we deliberately limit our attention to the n > 1 region 
of the phase diagram, as this is a typical situation in the 
experiments. 

We focus our attention on the behavior of the central 
compressibility, 



9/i 



(24) 



corresponding to the response of the density in the cen- 
tral region C to a variation of the chemical potential, 
driven exclusively by a change in the trapping potential 
(which leaves the parameters J and U unchanged) . Upon 
choosing a sufficiently narrow central region C satisfying 
the condition of Eq. [7j the central density is supposed to 
reproduce closely the behavior of the compressibility in 
the bulk system. Fig. [16] shows indeed that this is the 
case. In fact we observe that the central compressibility 
jumps from from a finite value to zero for critical values of 



the chemical potential which reproduce closely the shape 
of the Mott-insulating lobes of the bulk phase diagram 
of the Id Bose-Hubbard model. (The Mott insulating 
boundaries for the bulk system have been determined 
in Ref. |40j via density-matrix renormalization group for 
the n = 1 and n = 2 lobes, while we have determined 
the 71 = 3 lobe via quantum Monte Carlo). Remarkably, 
upon doubling the central region C from 10 to 20 sites 
(namely halving the resolution of the imaging in the ex- 
periment), the main features of the reconstructed phase 
diagram are only slightly altered, at least for what con- 
cerns the lobes n — 1, 2. This shows the robustness of 
the proposed approach under realistic experimental con- 
ditions. An interesting phenomenon to be observed in 



Fig. 16 is the "blue shift" (towards higher chemical po- 
tentials) of the n — 2 and n — 3 Mott lobes predicted 
by the vanishing of with respect to what is observed 
in a bulk system. This shift can be easily understood as 
a result of the finite width of the imaging region and as 
an effect of the finiteness of the central region exhibiting 
bulk physics, and, as we will discuss, it can be corrected 
for systematically. 

Blue shift of the lower lobe boundaries. The lower lobe 
boundaries are associated with a SF-MI transition in- 
duced by an increase of the chemical potential above 
a lower critical value iJ,c\j/U). A chemical potential 
fi > forces the atom density to an integer value, for 
which a Mott gap opens if J/U < {J/U)c corresponding 
to the lobe tip. When this phenomenon occurs in the 
center of the trap, the MI region which appears at first is 
extremely narrow in space (see Fig. |17[ a) for a sketch). 
Indeed it extends over a region of width Ar such that 
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FIG. 17: Mechanisms of blue shift of the Mott lobes obtained 
from the central compressibility, a) At the lower critical chem- 
ical potential fj,c'' a Mott region only forms on the central 
site, and one needs an increase of the chemical potential by a 
term to establish a Mott phase throughout the imaging 
region (pink-shaded), b) At the upper critical potential /ii"' 
of the bulk system, the trapped system has not yet devel- 
oped a superfluid core, because of the additional zero-point 
energy ezp{n) due to trapping in the harmonic oscillator po- 
tential. Only when the chemical potential reaches the blue- 
shifted value /ii"' -I- ezp('^), the first particle of the superfluid 
core appears (marked in red). 



the chemical potential modulation due to the trap over 
this region is lower than the MI hole gap Ah, namely 
Vt(Ar/2)^ « Ah- The central compressibility detects 
the appearance of an incompressible MI phase in the trap 
center only when Ar reaches the size |C| of the imaging 
region, namely for a critical bulk gap sa Vt.c{\C\/2y . 
Here Vt^c is the trapping potential corresponding to the 
effective chemical potential fi{Vt^c^ -^j J, U) at which the 
hole gap takes the value A^. This chemical potential is 
simply given by /xi' + AJ^. This implies that, for a strong 
confinement Vj, as that required to obtain n = 2 or n = 3 
MI lobes in the system, one observes the onset of a MI 
phase for a blue-shifted critical chemical potential 
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Remarkably this shift depends exclusively on experimen- 
tally known parameters Vt and |C|, so that one can cor- 
rect a posteriori for it when determining the phase dia- 
gram experimentally. 

Blue shift of the upper lobe boundaries. The upper 
lobe boundaries are associated with a MI-SF transition, 
induced by an increase of the chemical potential. Reach- 
ing an upper critical value the chemical potential 
closes the Mott gap of the MI with n particles per site. 



and it therefore adds the first particles on top of the 
Mott insulator. In a bulk system these particles are 
essentially added in a fc = state, and they undergo 
condensation. In a trap, the corresponding phenomenon 
would be the appearance, induced by trap squeezing, 
of the first particle over the central Mott plateau with 
n particles per site (see Fig. [iTjb)). At variance with 
the bulk case, these particles are not fully delocalized, 
but they still feel the overall trapping potential - hence 
they have a residual zero-point energy, due to the con- 
finement in the trap. To estimate this shift in the zero 
point energy, we can proceed in the spirit of the ex- 
tended fermionization approach |5S], (see Sec. 



IV and 



Appendix |Aj and use the approximation of representing 
the extra particle added in the trap center as moving 
on an inert background at fixed density n, whose only 
effect is to increase the hopping amplitude of the par- 
ticle from J to J(n -f 1). Hence the zero point energy 
can be estimated from the solution of a harmonic oscilla- 
tor of effective mass m* = /[2J{n + 1)] and frequency 
fujj = -\/4(n -f l)Vt J. This estimate is valid as long as the 
width of the central MI plateau is larger than the width 
of the ground state of the harmonic oscillator, namely 
a = [J {n-\-l) /VtY^ . Hence the shift in zero-point energy 
due to confinement is £zd (n) ^ huj/2 = ^{n + l)VtJ. 
This leads to a shift in the critical chemical potential 
needed to "dope" a particle onto the central MI plateau 
with respect to that needed in a bulk system; as a con- 
sequence, the upper critical chemical potential jic^ mea- 
sured by the central compressibility is blue-shifted as 
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Similarly to what observed for the lower lobe boundary, 
one can correct a posteriori for this shift using the ex- 
perimental parameters, and hence reconstruct the bulk 
value of /Xc"' . 

Fig. [18] shows the lower and upper critical chemical po- 
tentials of the n — 2 MI lobe, as extracted from the van- 
ishing of the central compressibility of iV = 100 bosons 
imaged over the central 20 sites, and compared with the 
numerically exact results of Ref. 40l The blue shift of 
the lower critical chemical potential is quite sizable, but 
it is seen to be very well corrected for by the prediction 
of Eq. ( 25 ) . The blue shift of the upper critical chemical 



potential is instead much weaker, in agreement with the 
prediction of Eq. (26). This shows that trap squeezing 



allows to reconstruct quite accurately the quantum phase 
transition lines of the Bose-Hubbard model, even in the 
case of a coarse imaging of a central region of ^ 20 sites 
in the trap center. 

Trap squeezing supplemented with imaging over |C| 
sites can reliably detect MI behavior only as long as the 
MI plateau occurring in the trap center can reach the 
same width as that of the central C region: this con- 
dition guarantees that (/ic'')' < (/ii"'')', namely that, 
upon increase the chemical potential, the (blue-shifted) 
lower chemical potential is reached before a new super- 
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FIG. 18: n = 2 Mott lobe of the Id Bose-Hubbard model 
as determined by trap squeezing of = 100 bosons, imaged 
over the central 20 sites. The open circles correspond to the 
bare values of the lower and upper critical chemical potentials 
extracted from the vanishing of the central compressibility for 
\C\ = 20 (Fig. 16 1; the closed circles corres pon d to t he values 
corrected for blue shift, according to Eqs. (|25[l and (26 1. The 



correction leads to a very good agreement with the DMRG 
results of Ref. gDl 



fluid core forms. The MI region can no longer be reliably 
observed if Vt reaches the value Vt satisfying the condi- 
tion (^i'V - ^"("^^ 



namely 



Vt (|C|/2)2 - ^{n+l)VtJ = - ^4^'> = Aph (27) 



where Aph is the particle-hole gap of the Mott insulating 
phase. The above condition gives 
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where R ~ \C\/2 ^ 1. The merging of the blue-shifted 
lower Mott lobe boundary with the upper one is seen to 
occur e.g. in Fig. [l8] for J/U > 0.12. To be able to 
observe the tip of the Mott lobe via trap squeezing, it 
would be necessary to increase the total particle number 
so that Vt can be decreased below Vt, maintaining the 
same chemical potential in accordance with Eq. (13). 



So far we have focused our discussion on the spectral 
properties associated with the compressibility. In the fol- 
lowing we will also discuss the evolution of the coherence 
properties under trap squeezing, as an alternative tool to 
investigate the phase diagram of the system. 



B. Time-of- flight observables and condensate 
compressibility 

At variance with diagonal observables such as density 
and compressibility, which can be measured locally, the 
phase coherence properties, encoded in the momentum 
distribution, are generally measured at the global level. 
We define the momentum distribution as 



n(fc) 
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(29) 



where the normalization to the boson number N is chosen 
so as to account for the actual spatial extent of the cloud 
|65) . The momentum distribution is probed via time of 
flight measurements, by collecting the interferogram of 
the matter wave emitted from all sites of the optical lat- 
tice. Given the strong inhomogeneity of the trapped sys- 
tem, one can suspect that the momentum distribution 
of a trapped gas is generally not a faithful portrait of 
the corresponding distribution of a homogeneous system 
with the same global chemical potential. Yet some re- 
marks are in order. For any given tolerance e, one can 
select a central region C satisfying Eq. ([7|, and a com- 
plementary region C. Hence the momentum distribution 
can be decomposed into three contributions, two of which 
involve the central region C 



i{k) = nc{k)+nc^c{k)+nc(k) (30) 



where 



and 
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We would like to argue in the following that the behav- 
ior of the global momentum distribution is strongly de- 
pendent on the contribution nc{k) from the central re- 
gion. Indeed the central region of the trap is the only 
simply connected region, which implies that the pairs of 
sites {ij) associated with this region are the ones at the 
shortest mutual distance \ri — rj\. This latter topological 
aspect has strong consequences on the coherence prop- 
erties of the global system. Indeed one can envision in 
general two situations: 1) If the central region C is lo- 
cally in a SF state, its coherence properties are going to 
give the dominant contribution to the sum of Eq. ( |29[ ), 
because nc{k) contains the largest number of strongly 
correlated pairs (ij) at short distance; moreover if the 
wings are SF too, the central region will also contribute 
with the tlq (j(k) term, which generally dominates over 
Ti(j{k) because the central region is at a higher density 
than its complement C. 2) If the central region is in a lo- 
cal MI phase, extending over a region of linear size ~ Rc, 
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-^/^ n{k = 0) -^/^ FWHM (a-^) 

FIG. 19: Evolution of time-of-flight observables for A*' = 100 bosons under variable trapping frequency. The left panel depicts 
the condensate number n^^^ , and the right panel the full width at half majcimum (FWHM) of the momentum distribution n^, 
in units of where a is the lattice spacing. 



there might still be a SF halo in the region C, but it will 
give a weak contribution to the n{k) because the pairs of 
sites (ij) involved in nQ{k) are separated by a distance 
~ Rc, and sites separated by the central Ml core are 
very weakly correlated, given that particles cannot co- 
herently propagate across C. An equivalent statement to 
the latter one is that in a d-dimensional system the SF 
halo surrounding a Ml region is effectively quasi- (d — 1)- 
dimensional, and consequently it has reduced coherence 
properties, due to the enhanced role of quantum fluc- 
tuations in lower dimensions. Hence it is reasonable to 
believe that the momentum distribution mainly reflects 
the coherence properties of the central C region, at least 
for one- and two-dimensional systems. 

Fig. [19] shows the height of the k = peak and the 
full width at half maximum (FWHM) of the n{k) in a 
system of = 100 bosons under trap squeezing con- 
trol of the chemical potential, and for various optical lat- 
tice depths. Plotting these quantities on the (J/J7, /x/C/) 
plane, one sees that sharp drops in the global coherence 
of the trapped system are observed close to the SF-Ml 
boundary lines of the bulk phase diagram of the Bose- 
Hubbard model. As in the case of the global compress- 
ibility, for high enough density the sharpest features in 
the coherence properties are observed for blue-shifted val- 
ues of the chemical potential [66]. This result seems to 
confirm our expectation that the global properties of the 
momentum distribution are strongly dominated by the 
behavior of the system at the trap center, which best 
approximates the bulk behavior. 

To make this connection even more quantitative, we 
can consider again the partition of the system into a 
central C region and a complementary peripheral C re- 
gion as a partition of the trapped gas into a system 
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FIG. 20; Absolute value of the condensate compressibility 
Kg = dn{k = 0)/dfi (in units of U~^) of N = 100 bosons 
under variable trapping frequency. 



and a particle reservoir, respectively (Fig. [T]). Trap 
squeezing increases the chemical potential of the reservoir 
with respect to that of the system, and it hence pumps 
particles from the former to the latter. If the system 
experiences (quasi-)condensation, a (quasi-)macroscopic 
fraction of the pumped particles will enter the (quasi- 
)condensate. Alternatively, if the central region has a 
density which approaches an integer value from below, 
nc < int{nc) + 1/2, condensation should be rather 
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thought of in terms of holes (of density int (nc) + 1 — nc); 
in that case a (quasi-)niacroscopic fraction of the holes 
removed from the center under trap squeezing belongs to 
the (quasi-)condensate. In both cases, if the trap center is 
in a local (quasi-)condensate state, one expects that the 
variations of the ri(fc = 0) peak in the momentum dis- 
tribution will be proportional to (for true condensation) 
or will go as a power law of (for quasi-condensation) the 
variation of the density in the trap center. On the other 
hand, if the center of the trap is in a MI insulating phase, 
this prevents the particle (hole) transfer, and it inhibits 
the variations of the momentum distribution. 

From the above arguments, one expects the behavior 
of the momentum distribution to be extremely sensitive 
to trap squeezing; this sensitivity can be quantified by 
the condensate compressibility 

Ko ^n{k = 0) (33) 

which can capture the (quasi-) macroscopic transfer of 
particles from the C reservoir into a condensed state in 
the center C (or viceversa for hole condensation). The 
condensate compressibility is positive for particle (quasi- 
) condensation in the trap center, and negative for hole 
(quasi-)condensation. Consequently |ko| is expected to 
be finite in presence of condensation in C (apart from 
the zeros corresponding to the passage from particle to 
hole condensation) , and to be essentially vanishing corre- 
sponding to MI behavior. Fig. [20] confirms this expecta- 
tions, showing that |ko| is vanishing corresponding to the 
bulk MI regions (apart from the usual blue shift of the 
critical chemical potentials) while it has generally finite 
values corresponding to the SF regions. The behavior in 
the SF regions is highly inhomogeneous, due to the fre- 
quent passage from particle to hole condensation. In any 
case, one can easily tell apart the accidental vanishing of 
the condensate compressibility in a SF regime from the 
systematic vanishing in the MI regime e.g. by inspection 
of the FWHM, which is much smaller than 27r/a in the 
SF case, and of the order of 27r/a in the MI case (a being 
the lattice spacing). 

VII. APPLICATIONS: BOSE GLASS IN A 
THREE-COLOR SUPERLATTICE 

A. Signatures of a Bose glass 

In the previous section we have seen that trap squeez- 
ing allows to faithfully reconstruct the phase diagram of 
the ordinary Bose-Hubbard model, either via the mea- 
surement of the diagonal properties in the central region 
of the trap, or via the measurement of the evolution of the 
global off-diagonal properties (momentum distribution). 
This latter aspect is due to the fact that spectral proper- 
ties and coherence properties can both characterize fully 
the two phases (SF and MI) present in the phase diagram 
of the Bose-Hubbard model. On the other hand, spectral 



and coherence properties become complementary pieces 
of information in presence of disorder in the system, in- 
ducing a novel phase in the phase diagram: the Bose 
glass ^8]. This extra phase is in fact characterized by a 
finite compressibility coexisting with short-range phase 
coherence, and it needs the joint measurement of local 
diagonal properties and global momentum distribution 
to be detected. 
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FIG. 21: Trap squeezing in the Bose glass phase. a) 
Schematic view of the density profile of a boson gas in a 
disordered optical lattice. 6^ In presence of a Bose glass be- 
havior in the trap center, upon infinitesimal trap squeezing 
Vt — >■ Vt+5Vt a particle can be transfered from the trap wings 
to a localized state in the center. Due to its localized nature, 
this additional particle does not contribute significantly to the 
coherent fraction, which remains essentially unchanged. 

In Ref. |30j we have shown that trap squeezing can re- 
veal the Bose glass phase appearing in the Bose-Hubbard 
model with a one-color quasi-periodic potential, real- 
ized by a two-color incommensurate optical superlattice 
[211 IMl EH [68]- The fingerprint of the Bose glass in a 
trap is the observation of a finite central compressibility 
Kc, coexisting with a weak condensate compressibility 
kq, and with the observation of a broad peak in the mo- 
mentum distribution. 

Fig. [21] sketches the response to trap squeezing of a 
Bose glass occurring in the trap center. By definition, 
the Bose glass phase admits particle excitations at ar- 
bitrarily low energy, which means that a small increase 
in the trap frequency Vt and, consequently, in the global 
chemical potential /x, is able to transfer a particle to the 
center of the trap, resulting in a finite central compress- 
ibility Kc- It is important to emphasize that this is only 
true on average over the disorder distribution: indeed, 
averaging over the disorder statistics, one expects always 
the existence of one or several disorder arrangements in 
the trap center which admit the transfer of a particle 
from the wings to the center at an infinitesimal energy 
cost. Explicit disorder averaging is crucial, because we 
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cannot expect in general the central region of the sys- 
tem to be self-averaging, due to its relatively small size. 
In the following we will consider quasi-periodic super- 
lattices, which admit random relative spatial phases be- 
tween the various Fourier components. In this case, the 
full disorder statistics in the trap is sampled by consid- 
ering different values of the relative phases. 

If the shift in chemical potential associated with trap 
squeezing is not bringing the system across a Bose-glass- 
to-superfluid transition, the particles that are transfered 
from the wings to the center will occupy localized quasi- 
particle states - meaning that the density variation that 
their transfer causes will be well localized in space. This 
implies that the new particle appearing in the center re- 
mains highly spread in momentum, and that it only con- 
tributes a small amount to the k = peak in the mo- 
mentum distribution. As a consequence, the coherence 
properties are only weakly altered by the particle trans- 
fer, and the condensate compressibility remains small. 

While a finite central compressibility kc is a solid cri- 
terion (and a pre-requisite) for the identification of a 
Bose glass, it is a priori more difficult to fix a criterion 
on the condensate compressibility which would distin- 
guish a Bose glass from a weakly correlated superfluid. 
In fact, in both cases the condensate compressibility is 
generally finite, given that a global increase of the den- 
sity leads to a global change of the momentum distribu- 
tion, regardless of whether the system exhibits condensa- 
tion or not. The only rigorous approach to discriminate 
a compressible insulator from a condensate is to study 
the scaling of the condensate occupation n{k = 0) un- 
der increasing particle number and decreasing trapping 
frequency, namely at fixed chemical potential; this will 
be the subject of a future publication. Nonetheless, a 
semi-quantitative criterion can be formulated based on 
the width of the momentum distribution. Indeed the in- 
verse of the FWHM of such a peak gives an estimate of 
the correlation length ^ in the system: in the following 
we will take ^ « 2 (FWHM) ~^ (based on the assumption 
that the momentum distribution n{k) can be approxi- 
mated by a Lorentzian with parameter ^). The system is 
obviously in an insulating state if ^ is far below the width 
of the atomic cloud. We can estimate such a width via 
the participation ratio 



PR 



(34) 



Hence the system can be safely considered as an insulator 
if e < PR. 



B. Three-color incommensurate superlattice 

We now focus our attention on the Id Bose-Hubbard 
in a two-color quasi-periodic potential, Eq. (17), real- 
ized via a three-color incommensurate superlattice. In 
Ref. [5S] we have shown that this potential, despite be- 



over a continuous region of the phase diagram. This is 
to be contrasted with the one-color quasi-periodic po- 
tential, realized via two-color superlattices [12], which 
features Bose glass behavior as well as incompressible 
incommensurate-band-insulator behavior [361 167] . ap- 
pearing in a tight alternation upon varying the chemical 
potential. In other words, while the one-color potential 
still retains features of its quasi-periodic nature, the two- 
color potential seems to mimic very closely the behavior 
of a purely random potential. 

In the following we focus indeed on the regime in which 
the system immersed in a two-color potential exhibits an 
extended Bose-glass phase, namely for potential strength 
V2 = U which removes completely the Mott insulator 
from the phase diagram. In particular we consider the 
case of a strong potential, V2 = 20J, which establishes 
a Bose glass for all the values of the chemical potential 
we have explored. The trapped version of the system 
is investigated in the case of = 100 bosons, and un- 
der averaging over the spatial phases (p, cj)' of the Fourier 
components of the potential, Eq. (17). In particular we 



investigate the behavior of the phase-averaged quanti- 
ties: the central density {nc)^,, the momentum distribu- 
tion {n{k))^ and its FWHM, and the participation ratio 
(PR)^. We will plot these quantities as a function of the 
phase-averaged chemical potential, which in Section HI 
has been shown to satisfy a simple scaling relationship 
as a function of the experimental parameters Vt and N . 
The data on the trapped system are compared with those 
of the bulk system, realized by a chain of L = 300 sites 
for a fixed value of the potential phases 0, 0'; in this 
case disorder averaging is provided by the extended size 
of the system, which shows to be fully self-averaging (we 
have explicitly checked this aspect by comparison with 
an even longer chain, L = 500). 



ing quasi-random, appears to support a Bose-glass phase 



Fig. 22 shows the evolution of the central density in 
the trapped system for a central region C of varying 
size, compared with the bulk behavior. The bulk be- 
havior shows a continuous variation of the density with 
the chemical potential, signaling a finite compressibility 
over the whole range we have considered. We observe 
that the phase-averaged central density closely follows 
the bulk behavior, and, in the case of TV = 100 bosons 
we consider, the dependence of the central density on the 
size of the C region is very weak up to densities n < 2, 
signaling that a faithful measurement of the n(/i) curve 
can be achieved at different imaging resolutions. To faith- 
fully reproduce the behavior at higher densities it is nec- 
essary to have a good resolution of around 10 sites, or, 
alternatively, to increase the number N of trapped parti- 
cles and reduce the trapping potential (thereby keeping 
the product N'^^^V^^^ constant) in order to increase the 
quasi-homogeneous region in the trap center. 

The insulating nature of the phase stabilized in the 
trap center is revealed by the momentum distribution. 
In Fig. [23] we show the evolution of the peak in the mo- 
mentum distribution under trap squeezing. The compar- 
ison with the data of the bulk system shows that the 
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FIG. 22: Phase-averaged central density of a system oi N = 
100 one-dimensional bosons trapped in a parabolic potential 



plus a three-color superlattice, Eq. (171 (see text for the pa- 
rameters). The results for the trapped system are compared 
with those of a bulk system with L = 300 sites. 



FIG. 24: Correlation length of system oi N = 100 one- 
dimensional bosons trapped in a parabolic potential plus a 
three-color superlattice (same as in Fig. |22| , estimated via the 
inverse FWHM of the momentum distribution. This quantity 
is compared with the effective cloud size, given by the phase- 
averaged participation ratio. 
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FIG. 23: fc = peak in the momentum distribution for a 
system of A'^ = 100 one-dimensional bosons trapped in a 
parabolic potential plus a three-color superlattice (same as 
in Fig. 22 1. The results are compared with those of a bulk 
system with L = 300 sites. 



easily drawn on the basis of the observation of the mo- 
mentum distribution. In fact, at variance with the Mott 
insulating state, the condensate compressibility kq ap- 
pears to be finite everywhere, although large excursions 
of the chemical potential {e.g. from 0.5 U to U) only 
slightly change the condensate population (by about 10% 
in the case mentioned). As already mentioned, the most 
dramatic signature of the insulating behavior would be 
revealed upon scaling the system size. Yet, in this case 
the insulating behavior is so pronounced that it is al- 
ready revealed at the level of the single system size we 
considered. In fact, as shown in Fig. [24j the correlation 
length estimated via the inverse of the full width at half 
maximum is considerably smaller than the system size, 
estimated via the participation ratio (it is more than an 
order of magnitude smaller for n/U < 2). Hence the ob- 
servation of the central density and its evolution under 
trap squeezing, together with the measurement of the 
momentum distribution, would experimentally allow to 
conclude unambiguously about the Bose-glass nature of 
the phase realized in this system. 



trapped system reproduces the main features of the bulk 
behavior, including some subtle ones (for instance, the 
anomaly for ^/U ~ 2); this further confirms that the 
momentum distribution is dominated by the behavior in 
the trap center, which best approximates the bulk be- 
havior. In the case of the bulk system, the inspection 
of the superfluid density [36] allows to conclude that the 
system is in an insulating state for all explored chemical 
potentials, as the superfluid density turns out to be iden- 
tically vanishing everywhere. This conclusion cannot be 



VIII. CENTRAL VS. GLOBAL 
COMPRESSIBILITY 

Recent experiments have probed for the first time the 
compressibility of a fermionic gas in an optical lattice 
via trap squeezing and in situ imaging of the cloud 
width [31]. The same technique is obviously applicable 
to bosons as well, and it has been proposed as a practi- 
cal tool for the detection of the Bose glass (SHj . Probing 
the response of the global width of the atomic cloud to 
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trap squeezing has the obvious advantage that the imag- 
ing resolution required in this case is much less than that 
demanded by the technique we propose. Here we would 
like to point out a few drawbacks, particularly associated 
with the detection of the Bose-glass phase. A discussion 
of the drawbacks for the case of repulsive fermions can 
be found in Ref. [5TJ 
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FIG. 25: Central density (averaged over 10 sites) of a system 
of = 100 trapped bosons, described by the one-dimensional 
Bose-Hubbard model with U = 20J in a parabolic trap. This 
quantity is compared with the effective global density of the 
atomic cloud, estimated via the inverse of the cloud radius R 
and via the inverse of the participation ratio. 



In Ref. [31] the cloud width is estimated via the cloud 
radius, namely the second moment of the density distri- 
bution (n,;) |70j (centered around the origin) 



R = 



(35) 



The response of the atomic cloud to the variation of the 
trap frequency is then quantified in terms of the variation 
in the cloud radius R; given that N/ R"^ gives an estimate 
of the effective average density of the cloud, one can de- 
fine a global cloud compressibility the derivative of 
with respect to the trapping potential 



1 dR 



(36) 



This is a measure of the compressibility of an inhomoge- 
neous trapped system, which collects contributions from 
all its various parts. In particular, no matter what state 
the system core is in, the dilute halo around the core is 
always in a compressible state, so that the global com- 
pressibility is necessarily always finite, although it will 
be significantly higher or smaller depending on whether 
the core is in a locally compressible state or not. Hence 
the realization of a specific bulk phase in the system core 



does not result in a qualitatively different behavior of the 
global compressibility, but only a quantitatively different 
one. This aspect makes the discrimination of the phases 
realized in the system core very hard, and even more so 
the location of the phase boundaries. In particular, as 
we will discuss below, the identification of some phases 
might even be completely missed. 



Fig. 25 shows the global density N/R for 100 bosons 
described by the Id Bose-Hubbard model with U = 20 J 
and for a variable trapping potential, controlling the 
chemical potential of the system. An alternative defini- 
tion of the global average density, which turns out to be 
much closer to the true average density, is iV/PR, based 
on the participation ratio PR, Eq. ( 34 1 , which gives the 



effective volume occupied by the system. The PR is very 
close to size of the region over which the atomic den- 
sity is non-zero, while the cloud radius R is closer to the 
FWHM of the density profile. Despite the apparent vari- 
ety in the definitions. Fig. 25 shows that N/R and A^/PR 



exhibit precisely the same behavior, and they essentially 
differ by a multiplicative factor. Hence we will use them 
interchangeably in the following. 



In Fig. 25 the global densities N/R and A^/PR are com- 
pared with the central density {nc) for a central region 
of 10 sites - as shown in Ref. t30j for the same parameter 
set, this latter quantity reproduces rather faithfully the 
density of the bulk system. For the chemical potential ex- 
cursion considered in this study, we observe three Mott 
plateaus, which are clearly manifested in the central den- 
sity, and which also show up in the global densities via 
a significant change of slope in the dependence on the 
chemical potential 71J. In this case one could argue that 
the global compressibility contains essentially the same 
information as the central one. At the same time, it is 
evident that important pieces of information are missing 
in the global compressibility: under the assumption that 
a loss in global compressibility is always equivalent to 
the appearance of an incompressible phase in the system 
core, it is quite hard to locate the phase boundaries be- 
tween core-compressible and core-incompressible phases, 
as the changes in slope are quite smooth (despite the 
commensurate-incommensurate transitions in the core 
being very sharp for what concerns the central compress- 
ibility, as we have already discussed in Sec. VI). More- 



over, the global densities contain very little information 
on the central density at which the core-incompressible 
phases are realized. As a general remark, one could say 
that the global density and global compressibility can 
give information on the alternation of phases and their 
nature provided that one has a prior knowledge of the 
topology of the phase diagram of the system, and of the 
particular phases that are realized there. 

A much harder (yet more general) situation is the one 
in which one does not know the succession of phases re- 
alized in the system - although the set of the possible 
phases might still be known. This is a generic situation 
in which one would like an analog quantum simulator 
to work. As an example of a complex phase diagram. 
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FIG. 26: Central density (averaged over 10 sites) of a system 
of A'^ = 100 bosons trapped in a parabolic potential plus a 
two-color incommensurate superlattice (see text). This quan- 
tity is averaged over different realizations of the relative spa- 
tial phase of the superlattice. Comparison is made with the 
global density, estimated via the phase-averaged inverse par- 
ticipation ratio. 



we consider the case of the Bose-Hubbard model in an 
incommensurate one-color potential [5^ 1^ [5H] . As al- 
ready mentioned elsewhere in this paper, this system ex- 
hibits four bosonic phases: Mott insulator, superfluid, 
Bose glass, and incommensurate band insulator. Follow- 
ing Ref. [33], we consider a system of = 100 bosons in 
a strong incommensurate potential V2 = U = 20J with 
incommensuration parameter a = 830/1076, and under 
variable trapping potential. Fig. [26] shows the phase- 
averaged global density of such a system (iV/PR)^, com- 
pared with the phase-averaged central density {nc)(t, for 
a region of 10 sites. As seen in Ref. 136 it his latter quantity 
reproduces most of the features of the density of a bulk 
system for low enough values, n < 2.5. In particular a 
distinctive feature is represented by plateaus at incom- 
mensurate densities, related to the potential parameter 
a, and corresponding to incommensurate band insulating 
phases. The plateau regions are separated by compress- 
ible regions manifesting either Bose glass or superfluid 
behavior (for such strong external potential the Mott in- 
sulator is completely washed out, and no plateau appears 
at integer densities). 

When inspecting the global density, one can still ob- 
serve a weak decrease of its slope as a function of the 
chemical potential corresponding to the incommensurate 
plateaus. Nonetheless, without prior knowledge of the 
bulk phase diagram, it is very hard to understand the 
nature of this less compressible behavior. First of all, 
it is a priori unclear whether a local minimum in the 
compressibility should always be associated with the ap- 
pearance of a strictly incompressible behavior in the core. 



and, even working under this assumption, the estimate of 
the phase boundaries is somewhat arbitrary. Moreover, 
it is not possible in principle to decide whether the less 
compressible phases correspond to ordinary Mott insula- 
tors or to some other form of incompressible insulating 
state, because one lacks definite information on the filling 
at the trap center. Hence we can conclude that the global 
density and compressibility do not allow in general to re- 
construct the phase diagram of the Bose-Hubbard model 
in a quasi-periodic potential. 



IX. CONCLUSIONS 

In this paper we have shown how the control on the 
trapping potential for bosons in optical lattices provides 
a very practical experimental knob on the total chemi- 
cal potential of the system. Indeed, even in the strongly 
correlated case of bosons described by the Bose-Hubbard 
model with or without an external superlattice potential, 
the global chemical potential at zero-temperature turns 
out to be related to the trapping potential (and to the 
other Hamiltonian parameters) via a simple scaling rela- 
tion. This relation takes the form of a modified Thomas- 
Fermi scaling relation, corrected by quantum effects, and 
valid in the case of large enough densities or weak enough 
interactions, in which lattice commensuration effects are 
negligible. Corrections beyond this simple expression are 
quantitatively captured at the level of a simple lattice 
mean- field theory (Gutzwiller Ansatz), which agrees sur- 
prinsingly well with quantum Monte Carlo results. The 
possibility of tuning the chemical potential via the trap 
gives access to the exploration of the phase diagram 
of Hamiltonian models in the grand-canonical ensemble, 
even if cold-atom experiments are generally performed at 
a fixed number of particles. 

If supplemented with an imaging technique able to 
measure selectively the average density nc in a finite re- 
gion C around the trap center, the a priori knowledge of 
the global chemical potential /i of the system in an op- 
tical lattice allows to reconstruct the density-ws.-/i curve 
of the model implemented in the experiment without the 
trap. In particular this gives access to the bulk com- 
pressibility At, estimated via the central compressibility 
He — One /d 11. We show that central compressibility 
measurements tell apart very clearly Mott insulator from 
superfluid phases, and hence they allow to fully recon- 
struct the phase diagram of the Bose-Hubbard model at 
zero temperature in the grand-canonical ensemble. The 
addition of an incommensurate superlattice |22j intro- 
duces Bose glass regions in the phase diagram, whose 
unambiguous detection can be achieved via the joint mea- 
surement of the central compressibility and of the global 
coherence of the system, as well as its response to trap 
squeezing. 

The broader perspective of a truly useful quantum sim- 
ulation with cold atoms requires to implement computa- 
tionally more challenging models than lattice bosons - 
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e.g., to simulate lattice fermions. To do so with minimal 
"classical input" , on needs a strategy to know in advance 
all the parameters of the system without the necessity 
an ab-initio classical simulation of the model of interest. 
While it is still unclear how to gain a priori knowledge on 
the temperature without extensive classical simulations 
pn [T5] , the case study of the Bose- Hubbard model, pre- 
sented in this work, suggests that mean-field variational 
Ansatzes might be successful to reconstruct a priori the 
chemical potential of strongly correlated system in their 
ground state, even when the same variational Ansatzes 
do not provide an accurate description of the microscopic 
behavior of the system. 
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Appendix A: Jordan- Wigner fermionization and 
extended fermionization 



Hardcore boson operators \bi, foj] = Q {i ^ j), {bi, b\} ^ 

1 can be mapped onto spinless fermions /i, // via the 
Jordan- Wigner transformation [43 . 



fe=i 



fe=i 



so that the Hamiltonian of infinitely repulsive bosons in 
an external potential Vi becomes that of non-interacting 
spinless fermions: 



(A2) 



for which we assume open boundary conditions. 

This approach has been recently extended to the case 
of strongly repulsive Id softcore bosons in a trap [55] for 
U/J > {U/J)c,n, where {U/J)c,n is the critical value for 
the superfiuid-to-Mott-insulator transition at the largest 
integer density n reproduced locally in the trap. Under 
this assumption, the density profile in the trap takes a 
characteristic "wedding-cake" structure, namely regions 
at densities in the interval [m — l,m] (with integer m) 
are separated from regions with densities in the interval 
[m, m + 1] by a sizable Mott plateau at integer density 
m. If U/J is sufficiently large, density fluctuations in 
the plateau can be neglected, which means that density 
fluctuations in the [m — l,m\ region are essentially de- 
coupled from those in the [to, to -|- 1] region, and that 
the lattice sites in each of these regions only experience 
density fluctuations between the two extremal values of 
the average density. This fundamental restriction of the 
local Hilbert space to a set of two-level systems implies 



Appendix B: Thomas- Fermi results for the atom 
number distribution in low-dimensional systems 

We consider a gas of weakly interacting bosons in a 
three-dimensional trap of frequencies (w^; , cj^ , ) and in 
a standing wave potential along 3 — d spatial dimen- 
sions, deflning an array of d-dimensional trapping ele- 
ments (tubes for d = 1, pancakes for d = 2). In the 
weakly interacting case, we describe the system within 
Gross-Pitaevskii (GP) theory in continuum space for d 
dimensions, and on a lattice for the remaining 3 — d ones. 
To calculate the density profile of the atoms, we make 
use of the Thomas- Fermi approximation, which amounts 
to minimizing the potential-energy part of the GP energy 
functional. 



1. d = 1 (tubes) 

In this case the system is immersed in two standing 
waves along the x and y directions, with wavelength A^; 
and Aj, respectively, defining a 2d array of tubes. Being 
dx — (da;,0), dy = {0,dy) the vectors of the 2d array, 
with da — Aq/2 (a = x,y), we identify the position of 
the tubes as R{ix, iy) = ix dx+iy dy. Imagining that the 
atoms can only occupy the lowest Wannier states of each 
tube in the x and y directions {wx and Wy, respectively), 
we can discretize the space in those directions, so that 
the GP wavefunction can be written as 5* — {ix , iy', z) . 
The potential energy part of the associated GP functional 
reads: 



E, 



(pot) 
GP 



-gild) 



'9{ix,iy;z)\'^ 



^xil + eyil + -mulz^ - n^^'^^ ] |«'(i„ z^; z)\' 



where — muy^d1^/2, and 
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g dx dy \wx{x)\ \wy{y)\ 



(Bl) 



(B2) 
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is the effective coupling constant for particles trapped in 
the lowest Wannier states in the x and y direction, related 
to the bare coupling constant 47rS^ a^/m (a^ — s-wave 
scattering length) . Taking a Gaussian approximation for 
the Wannier states 1721 one obtains 



array of pancakes. In analogy with what done in the 
previous section, we define a discretized GP wavefunc- 
tion = ^{x,y]iz)^ whose associated potential energy 
is described by the GP functional 



-(irf) _ kxky f VoxVpy 



1/4 



(B3) 



4p°"[*,v1/*]=^ dxdy ^-—\^{x,y,i,t + 



where ka = 2'K/Xa, Voa is the lattice depth in the a 
direction, and Era = fi^k^/{2m) the associated recoil 
energy. 

The Thomas- Fermi approximation gives for the num- 
ber of atoms in each tube the expression 



(Id) 



(B4) 

where the chemical potential [i^^'^^ is consistently given 
by: 



where we have introduced the effective coupling constant 



1/4 



(B8) 
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2 \ 3/2' 
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with u) = (w^cjjjWz)^/^. Integrating in Eq. (B4) on the re- 
gion in which the integrand is positive, one finally obtains 
the expression [75] : 



4^/2 



/Id) 



3/2 



(B6) 



with evident meaning of the symbols. 

Integrating the Thomas-Fermi expression for the den- 
sity within each pancake, we obtain the atom number 
distribution 



with chemical potential 



2d) _ , ,-2 



(B9) 



2. d — 2 (pancakes) 

In this case the system is immersed in a standing wave 
along the z direction with wavelength X^, defining a Id 



15 3(2d)^^^ /_,-,2x3/2' 



Stt 



-1 2/5 



(BIO) 



[1] I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 

80, 885 (2008). 
[2] M. Lewenstein, A. Sanpera, V. Ahufinger, B. Damski, A. 

Sen (De), and U. Sen, Adv. Phys. 56, 243 (2007). 
[3] R. P. Feynman, Int. J. Theor. Phys. 21, 467 (1982). 
[4] G. Murthy, D. P. Arovas, and A. Auerbach, Phys. Rev. 

B 55, 3104 (1997). 
[5] M. Boninsegni, J. Low Temp. Phys. 132, 39 (2003). 
[6] P. Sengupta, L. P. Pryadko, F. Alet, M. Troyer, and G. 

Schmid, Phys. Rev. Lett. 94, 207202 (2005). 
[7] V. W. Scarola, E. Demler, and S. Das Sarma, Phys. Rev. 

A 73, 051601 (2006). 
[8] M. P. A. Fisher, P. B. Weichman, G. Grinstein, and D. 

S. Fisher, Phys. Rev. B 40, 546 (1989). 
[9] R. B. Diener, Q. Zhou, H. Zhai, and T.-L. Ho, Phys. Rev. 

Lett. 98, 180404 (2007). 
[10] T.-L. Ho and Q. Zhou, Phys. Rev. Lett. 99, 120404 

(2007). 

[11] B. Capogrosso-Sansone, E. Kozik, N. Prokofev, and B. 



Svistunov, Phys. Rev. A 75, 013619 (2007). 
[12] L. PoUet, C. Kollath, K. Van Houcke, and M. Troyer, 

New J. Phys. 10 065001 (2008). 
[13] D. M. Weld, P. Medley H. Miyake, D. Hucul, D. E. 

Pritchard, and W. Ketterle, Phys. Rev. Lett. 103, 245301 

(2009). 

[14] D. McCay M. White, and B. DeMarco, Phys. Rev. A 79, 
063605 (2009). 

[15] D. McKay B. DeMarco, arXiv:0911.4143 (2009). 

[16] Q. Zhou and T.-L. Ho, arXiv;0908.3015 

[17] S. Trotzky, L. Pollet, F. Gerbier, U. Schnorrberger, I. 
Bloch, N.V. Prokofev, B. Svistunov, and M. Troyer, 
|arXiv:0905.4882| (2009). 

[18] R. Joerdens, L. Tarruell, D. Greif, T. Uehlinger, N. 
Strohmaier, H. Moritz, T. EssHnger, L. De Leo, C. Kol- 
lath, A. Georges, V. Scarola, L. Pollet, E. Burovski, E. 
Kozik, and M. Troyer, arXiv;0912.3790 (2009). 

[19] Q. Zhou, Y. Kato, N. Kawashima, and N. Trivedi, Phys. 
Rev. Lett. 103, 085701 (2009). 



27 



[20] T.-L. Ho and Q. Zhou, Nature Phys. 6, 131 (2009). 
[21] T. Stoferle, H. Moritz, C. Schori, M. Kohl, and T. 

Esslinger, Phys. Rev. Lett. 92, 130403 (2004). 
[22] L. Fallani , J. E. Lye, V. Guarrera, C. Fort, and M. In- 

guscio, Phys. Rev. Lett. 98, 130404 (2007). 
[23] D. Clement, N. Fabbri, L. Fallani, C. Fort and M. Ingus- 

cio, Phys. Rev. Lett. 102, 155301 (2009). 
[24] P. T. Ernst, S. Gotze, J. S. Krauser, K. Pyka, D.-S. 

Liihmann, D. Pfannkuche, K. Sengstock, Nature Phys. 

6, 56 (2009). 

[25] K. D. Nelson, X. Li and D. S. Weiss, Nature Phys. 3, 556 
(2007). 

[26] T. Gericke, P. Wiirtz, D. Reitz, T. Langen, and H. Ott, 
Nature Phys. 4, 949 (2008). 

[27] N. Gemelke, X. Zhang, C.-L. Hung, and C. Chin, Nature 
460, 995 (2009). 

[28] W. S. Bakr, J. I. Gillen, A. Peng, S. Foiling, and M. 
Greiner, Nature 462, 7269 (2009). 

[29] I mmanuel Bloch's group website: http: / / www.quantum- 
m unich.de /research / single-site- addressing- in-optical- 
lattices/ 

[30] T. Roscilde, New J. Phys. 11, 023019 (2009). 

[31] U. Schneider, L. Hackermuller, S. WiU, Th. Best, I. 

Bloch, T. A. Costi, R. W. Helmes, D. Rasch, and A. 

Rosch, Science 322, 1520 (2008). 
[32] S. Foiling, S. Trotzky, P. Cheinet, M. Feld, R. Saers, 

A. Widera, T. Miiller, and L Bloch, Nature 448, 1029 
(2007). 

[33] M. Anderlini, P.J. Lee, B.L. Brown, J. Sebby-Strabley, 
W.D. Phillips and J. V. Porto, et al, Nature 448, 452 
(2007). 

[34] S. Bergkvist, P. Henelius, and A. Rosengren, Phys. Rev. 

A 70, 053601 (2004). 
[35] S. Wessel, F. Alet, M. Troyer, and G. G. Batrouni, Phys. 

Rev. A 70, 053615 (2004). 
[36] T. Roscilde, Phys. Rev. A 77, 063605 (2008). 
[37] Y.R.P. Sortais, H. Marion, C. Tuchendler, A.M. Lance, 

M. Lamare, P. Fournet, C. Armellin, R. Mercier, G. 

Messin, A. Browaeys, and P. Grangier, Phys. Rev. A 75, 

013406 (2007). 

[38] C. J. Pethick and H. Smith, Bose-Einstem Condensation 

in Dilute Gases, Cambridge, 2002. 
[39] A. W. Sandvik, Phys. Rev. B 59, R14157 (1999); O. F. 

Syljuasen, Phys. Rev. E 67, 046701 (2003). 
[40] M. Rizzi, D. Rossini, G. De Chiara, S. Montangero, and 

R. Fazio, Phys. Rev. Lett. 95, 240404 (2005). 
[41] M. Rigol, A. Muramatsu, G. G. Batrouni, and R. T. 

Scalettar, Phys. Rev. Lett. 91, 130403 (2003). 
[42] M. Rigol and A. Muramatsu, Phys. Rev. A 69, 053612 

(2004); Opt. Commun. 243, 33 (2004). 
[43] E. Lieb, Schuhz, and D. Mattis, Ann. Phys. (NY) 16, 

406 (1961). 

[44] L. Fallani and C. Fort, private communication. 

[45] B. Capogrosso-Sansone, S. G. Syler, N. Prokof'ev, and 

B. Svistunov, Phys. Rev. A 77, 015602 (2008). 

[46] B. Capogrosso-Sansone, N. V. Prokof'ev, and B. V. Svis- 
tunov, Phys. Rev. B 75, 134302 (2007). 

[47] We notice that the fj, data for the d = 3 system in a su- 
perlattice are much more scattered than for all the other 
cases, as shown on the right panel of Fig. [7] This is due 
to strong size effects on the system considered, which are 
inevitable in d = 3 due to the short linear size of each 
spatial dimension (we have considered cubes with up to 
L = 18 sites on each side, but the particles are actu- 



ally confined on a region with linear extent which is even 
smaller). In presence of a superlattice potential, a system 
on such a finite size experiences only a few wells of the 
potential in each spatial dimension, and hence the nonlin- 
earities due to such potential become relevant and cannot 
be simply averaged out as in the TF/AL approach. 
[48] G. G. Batrouni, H. R. Krishnamurthy, K. W. Mahmud, 
V. G. Rousseau, and R. T. Scalettar, Phys. Rev. A 78, 
023627 (2008). 

[49] M. Rigol, G. G. Batrouni, V. G. Rousseau, and R. T. 

Scalettar, Phys. Rev. A 79, 053605 (2009). 
[50] R. Jordens, N. Strohmaier, K. Gunter, H. Moritz, and T. 

Esslinger, Nature 455, 204 (2008). 
[51] V. W. Scarola, L. PoUet, J. Oitmaa, and M. Troyer, Phys. 

Rev. Lett. 102, 135302 (2009). 
[52] T. Kinoshita, T. Wenger, and D. S. Weiss, Nature 440, 

900 (2006). 

[53] M. Greiner, O. Mandel, T. W. Hansch, and I. Bloch, 

Nature 419, 51 (2002). 
[54] J. Dziarmaga, Phys. Rev. Lett. 95, 245701 (2005). 
[55] G. PupiUo , A. M. Rey , C. J. Wilhams, and C. W. Clark, 

New J. Phys. 8, 161 (2006). 
[56] We estimate gaps as small as ~ 10"* J, and whose small- 

ness might be limited by the resolution of the grid with 

which we sample the trapping strength. 
[57] A. Aubry and C. Andre, Proc. Isr. Phys. Soc, ed. C.G. 

Kuper 3, Adam Hilger, Bristol, 1979. 
[58] J. B. Sokoloff, Phys. Rep. 126, 189 (1985). 
[59] L B. Spielman, W. D. Phillips, and J. V. Porto Phys. 

Rev. Lett. 98, 080404 (2007). 
[60] D. J. Han, M. T. DePue, and D. S. Weiss, Phys. Rev. A 

63, 023405 (2001). 
[61] W. Ketterle, K. B. Davis, M. A. Joffe, A. Martin, and D. 

E. Pritchard, Phys. Rev. Lett. 70, 2253 (1993). 
[62] C.-S. Chuu, F. Schreck, T. P. Meyrath, J. L. Hanssen, G. 

N. Price, and M. G. Raizen, Phys. Rev. Lett. 95, 260403 

(2005). 

[63] Y. Shin, C.H. Schunck, A. Schirotzek, and W. Ketterle, 
Nature 451, 689 (2008). 

[64] Y.-a. Liao, A. S. C. Rittner, T. Paprotta, W. Li, G. B. 
Partridge, R. G. Hulet, S. K. Baur, and E. J. Mueller, 
arXiv:0912.0092 

[65] The standard definition of the lattice momentum distri- 
bution would involve a normalization to the volume L'* 
of the quantization box, but in the case of a trapped sys- 
tem the quantization box is arbitrary. As an alternative, 
here we imagine that the quantization box is fixed to the 
number of bosons TV - which is a convenient convention, 
given that this number is kept fixed in our study of a 
trapped system. In this case the condensate occupation 
n{k = 0) takes the meaning of the number of particles 
with momenta in the box \k\ < ir/{aN), where a is the 
lattice spacing. 

[66] Ref. |49] has documented that local density signatures of 
the onset of the Mott phase in a trap show up at a higher 
U/ J value than in the bulk for one- and two-dimensional 
systems, and that this blue shift is much more pro- 
nounced for higher trapping potentials Vj. A similar 
shift has been experimentally documented in K. Jimenez- 
Garci'a, R. L. Compton, Y.-J. Lin, W. D. Phillips, J. V. 
Porto, and I. B. Spielman, arXiv: 1003. 1541 , based on 
time-of-fiight measurements of bosons trapped in quasi- 
two-dimensional optical lattices. We find that these ob- 
servations are fully consistent with our observations for 



28 



the central compressibility (Fig. 16 1 and for the time-of- 
flight observables (Fig. |19[ ) of a one-dimensional system, 
and that they can be easily explained with the fact that 
the onset of a Mott insulating phase with a small gap A is 
hindered by the trap inhomogeneity: indeed the trapping 
potential discards the Mott insulating behavior beyond 
a length scale \ri — ro| from the trap center, such that 
K(r. -ro)"~ A. 
[67] G. Roux, T. Barthel, I. P. McCuUoch, C. KoUath, U. 
SchoUwck, and T. Giamarchi, Phys. Rev. A 78, 023628 
(2008). 

[68] X. Deng, R. Citro, A. Minguzzi, and E. Orignac, Phys. 

Rev. A 78, 013625 (2008). 
[69] D. Delande and J. Zakrzewski, Phys. Rev. Lett. 102, 

085301 (2009). 

[70] In the case of Ref. j3Tj {rii) is actually the column- 
integrated density of a three-dimensional system, but this 



aspect is actually irrelevant for our discussion). 

[71] We avoid extracting the central and global compressibili- 
ties by differentiating the related densities, because of the 
significant numerical noise associated with the derivative. 
Nonetheless the changes in slope in Fig. [25] are evident. 

[72] W. Zwerger, Journ. of Optics B: Quantum Semiclass. 
Opt. 5, S9 (2003). 

[73] A similar formula appears in B. Paredes, A. Widera, V. 
Murg, O. Mandel, S. Foiling, I. Cirac, G. V. Shlyap- 
nikov, T. W. Hansch, and I. Bloch, Nature 429, 277 
(2004). That formula matches with ours for the case 
u]^ = u)y = u]^ if the factor of 57V/(27r7V(0, 0)) appear- 
ing there is changed into its inverse, 27rA'^(0, 0)/(57V) . We 
believe our formula to be correct, as it satisfies the sum 
rule i N^-^'^\ix,iy) ~ N, which is not satisfied in 
turn by the formula of Paredes et al. 



